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Abstract 


This  is  the  final  report  detailing  the  seed  project  research  which  was  funded  by  the 
U.S.  Army  through  its  European  Research  Office  during  the  period  September  1997  — 
November  1998. 

The  primary  purpose  of  the  research  funding  was  to  enable  Prof.  J.R.  Whiteman 
and  Dr.  S.  Shaw  (BICOM,  Brunei  University,.  Uxbridge,  England)  to  collaborate  with 
Dr.  A.R.  Johnson  (Army  Research  Laboratory,  Vehicle  Technology  Center,  NASA,  Lang¬ 
ley,  VA,  USA)  toward  developing  a  framework  for  the  adaptive  finite  element  solution  of 
quasistatic  viscoelasticity  problems. 

The  major  goal  was  to  develop  this  framework  in  the  context  of: 

•  the  practical  utility  of  the  internal  variable  formulation,  as  used  by  Dr.  Johnson; 

and 

•  the  theoretical  utility  of  the  hereditary  integral  formulation,  as  used  at  BICOM. 

The  first  of  these  allows  for  practical  software  to  be  developed  in  a  straightforward  way 
from  existing  linear  elasticity  codes,  while  the  second  facilitates  the  derivation  of  mathe¬ 
matically  rigorous  a  posteriori  error  estimates — the  essential  building  block  for  adaptive 
finite  element  solvers. 

During  the  project  we  proposed  and  developed  a  hybrid  algorithm  blending  the  best 
features  of  these  two  approaches.  Also,  we  implemented  our  a  posteriori  error  estimates  to 
produce  software  capable  of  automatic  spatial  error  control  based  on  adaptive  meshing.  A 
prototype  version  of  this  software  is  now  mounted  on  Dr.  Johnson’s  workstation  at  NASA, 
Langley,  and  a  short  “manual”  illustrating  the  configuration  and  use  of  this  software  is 
included  later  in  this  report.  Pull  details  of  the  work  undertaken  on  the  Seed  Project  are 
also  contained  in  the  following  pages. 
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Chapter  1 


Introduction 


1.1  The  report  in  a  nutshell 


This  report  is  based  on  a  research  logbook  (initiated  on  October  4,  1997)  of  our  collabo¬ 
rative  research  with  Dr.  Arthur  R.  Johnson  (at  ARL,  VTC,  NASA,  Langley)  on  Internal 
Variable  methods  applied  to: 

•  structural  viscoelasticity  problems  of  quasistatic  type  with  a  view  to  their  formula¬ 
tion  and  subsequent  adaptive  solution. 

The  research  was  sponsored  by  the  United  States  Army  through  a  European  Research 
Office  Seed  Project  (purchase  order  N68171-97-M-5763)  and  we  would  like  to  take  this 
opportunity  to  gratefully  acknowledge  this  support. 

This  first  part  of  the  document  contains  a  short  introduction  to  the  constitutive  mod¬ 
els  used  in  viscoelasticity  theory  in  terms  of  hereditary  integrals  (Chapter  2)  as  well  as 
internal  variables  (Chapters  2  and  3).  We  then  indicate  (Chapter  4)  how  the  hereditary 
integrals  lead  to  memory  terms  appearing  in  the  governing  equations,  thus  resulting  in 
Partial  Differential  Volterra  (or  pdv)  equation  problems  that  need  to  solved  in  order 
to  model  the  physics.  This  material  is  based  on  the  paper  [50].  We  also  include  a  small 
amount  of  detail  on  the  phenomena  of  Non-Fickian  Polymer  diffusion  and  some  recent 
attempts  by  applied  mathematicians  to  construct  mathematical  models.  In  parallel  with 
this  Seed  Project  we  at  BICOM  have  also  been  actively  engaged  in  developing  prototype 
software  to  model  these  effects.* 

The  Volterra  operator  (i.e.  the  memory)  can  be  removed  from  the  differential  equa¬ 
tions  by  appealing  to  the  viscoelastic  internal  variables,  in  the  manner  of  Johnson  and 
Tessler  in  [24]  for  example,  and  so  we  explore  also  (Chapter  5)  the  equations  resulting 
from  this  approach. 

By  “mixing  and  matching”  the  various  forms  of  the  constitutive  relationship  with  the 
governing  equations,  we  illustrate  in  the  next  part  of  the  report  that  it  is  possible  to  derive 
a  variety  of  differential  equation  problems  which  can  then  be  addressed  numerically  in 

*  The  computational  modelling  of  problems  of  nonlinear  diffusion  in  the  contexts  of  controlled  drug  release 
technology  and  percutaneous  drug  absorbtion  (Project  QR97/2,  ongoing) 
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order  to  model  the  problems  at  hand  (Chapter  6).  This  is  important  because,  as  we  also 
illustrate,  the  numerical  discretizations  of  various  equivalent  continuous  problems  can 
result  in  differing  numerical  solutions.  Prior  to  this  Seed  Project  Shaw  and  Whiteman 
had  already  spent  a  great  deal  of  effort  in  deriving  a  posteriori  error  estimates  for  the 
Volterra  formulation  of  the  quasistatic  problem  (see  [47,  48,  52,  57,  54,  55]),  and  this 
presents  a  natural  entry  point  into  the  task  of  providing  adaptive  software  based  on  error 
bounds  which  are  computable  in  terms  of  the  discrete  solution.  On  the  other  hand, 
Johnson  and  colleagues  have  expended  similar  amounts  of  effort  on  practical  solution 
software  based  on  internal  variable  formulations  [25,  24,  23,  26].  Since  these  will  result  in 
different  numerical  solutions  the  a  posteriori  estimates  for  the  Volterra  formulation  cannot 
be  used. 

Once  the  implications  of  this  had  been  fully  appreciated  (i.e.  a  duplication  of  effort 
in  order  to  derive  a  posteriori  estimates  for  internal  variable  formulations),  one  of  our 
foremost  goals  therefore  came  to  be  to  derive  a  “hybrid”  formulation  of  the  problem  that 
would  have  the  same  numerical  solution  as  the  discrete  Volterra  equations,  but  could  be 
implemented  algorithmically  in  the  same  manner  as  an  internal  variable  method.  Our 
proposed  algorithm  is  detailed  in  Chapter  7,  along  with  some  numerical  tests  to  illustrate 
that  the  numerical  scheme  does  indeed  behave  as  it  should  (in  terms  of  convergence  to  the 
continuous  solution  as  the  discretization  is  indefinitely  refined). 

We  move  on  to  adaptivity  in  Chapter  8.  We  first  recapitulate  the  a  posteriori  error 
estimate  as  given  in  [55]  and  then  use  it  to  implement  adaptive  meshing  in  the  linear  elas¬ 
ticity  context.  Once  the  technique  has  been  detailed  the  extension  to  the  time  dependent 
linear  quasistatic  viscoelasticity  case  is  straightforward,  and  we  include  several  examples 
to  demonstrate  adaptive  space-discretization-error  control  by  selective  mesh  refinement. 

As  noted  in  [47]  (and  seemingly  also  implied  by  others  in,  for  example,  [28,  2])  the 
robust  control  of  the  time  discretization  error  by  varying  the  time  steps  is  a  difficult 
issue  for  this  quasistatic  problem,  and  this  is  due  to  there  being  no  time  derivative  in 
the  governing  equation.  As  a  result  we  feel  that  guaranteed  temporal  error  control  can  be 
achieved  only  in  a  more  abstract  way  by  measuring  the  error  in  a  negative  norm.  Details  of 
our  work  in  this  direction  are  in  [53,  56],  the  first  of  which  gives  numerical  results  for  a  very 
simple  case  which,  nevertheless,  demonstrates  that  the  resulting  error  control  is  effective. 
The  extension  of  these  results  to  the  space-time  problem  introduced  and  considered  below 
is  non-trivial  and,  although  we  are  progressing,  no  results  are  yet  available. 

In  Chapter  9  we  give  a  brief  “manual”  describing  how  to  obtain,  configure  and  run  the 
software  that  was  used  to  produce  the  numerical  results  in  this  report.  The  report  closes 
with  Chapter  10,  In  which  we  suggest  further  research  work  that  would  follow  naturally 
from  this  Seed  Project,  and  then  we  give  a  complete  bibliography. 

The  remainder  of  this  chapter  is  devoted  first  to  a  brief  discussion  of  the  generic  PDV 
equations  that  arise  when  modelling  viscoelastic  effects,  and  then  to  establishing  a  context 
and  the  basic  notation  for  the  entire  report. 
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1.2  The  context 

In  their  simplest  form  viscoelastic  materials  exhibit  behaviour  characteristic  of  both  clas¬ 
sical  Hookean  solids  and  Newtonian  fluids.  The  resulting  effects  are  important  when  the 
material  is  deforming  under  an  applied  load.  This  load  could,  for  example,  be  due  to  exter¬ 
nally  applied  forces;  internal  deformation  caused  by  a  diffusing  penetrant;  or,  constrained 
thermal  expansion  caused  by  temperature  gradients.  See  for  example  [30,  7,  38].  Moreover, 
the  material  somehow  keeps  a  record  of  its  response  history  and,  for  this  reason,  viscoelas¬ 
tic  materials  axe  said  to  possess  memory.  This  memory  is  manifest  in  the  constitutive 
relationship  between  the  stress  and  strain  tensors,  a  and  £,  and  as  a  result  mathemat¬ 
ical  models  of  viscoelastic  behaviour  take  the  form  of  partial  differential  Volterra  (pdv) 
equation  problems.  The  canonical  forms  of  these  equations  axe:  the  elliptic  Volterra 
problem, 


Au(t)=f(t)+  [  B(t,s)u{s)ds ; 
Jo 


the  parabolic  Volterra  problem, 


(1.1) 


u'(t)  +  Au(t)  =  /(f)  +  [  B(t,  s)u(s)  ds ;  (1.2) 

Jo 

and,  the  hyperbolic  Volterra  problem, 

u"{t)  +  Au(t)  =  f{t)  +  f  B(t ,  s)u(s)  ds.  (1.3) 

Jo 

These  are  supplied  with  initial  and/or  boundary  data  as  appropriate,  and  the  dependence 
on  the  space  variable  x  is  suppressed.  In  these  problems  we  use  A  and  B(t,  s)  to  represent 
partial  differential  operators  (acting  only  in  the  space  variables)  where,  for  example,  we 
could  have 

A  :=  — V2  and  B(t,  s)  :=  — V  •  4>(ty  s)V, 

although  for  (1.1)  and  (1.3)  the  appropriate  form  for  A  is  the  linear  elasticity  operator — 
with  B(t,s)  “similar”. 

The  purpose  of  the  first  chapters  is  to  illustrate  how  the  memory  terms  arise  in  these 
equations  and  also  to  summarize  the  various  PDV  equations  used  when  modelling  problems 
of  quasistatic  and  dynamic  viscoelasticity,  and  non-Fickian  diffusion  in  polymers.  We 
also  indicate  some  of  the  numerical  analysis  work  that  has  been  carried  out  for  these 
problems  (but  we  do  not  claim  to  be  exhaustive,  for  a  fuller  account  see  [49]). 

Throughout,  the  positive  real  number  T  will  denote  a  final  time  and  we  use  J  :=  [0,  T ] 
and  X (0,T]  to  denote  time  intervals.  Also,  for  n  =  1,2  or  3  we  consider  f2  C  Rn  to  be 
an  open  bounded  domain  with  boundary  dCl.  Furthermore,  we  consider  dfl  in  the  form 


<9f2  :=  To  u  rv  with  r£>nrv  =  0, 
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where  the  closed  set  IV  C  5ft  is  called  the  Dirichlet  boundary  and  is  of  positive  measure 
so  that 

[  dr>o. 

JrD 


We  call  the  (possibly  empty)  open  set  IV  C  5f 2  the  Neumann  boundary.  The  reason 
for  this  terminology  is  the  obvious  one  where  we  refer  to  the  type  of  boundary  condition 
specified  on  these  subsets.  We  indicate  vector-valued  quantities  with  boldface  so  that,  for 
example,  we  use  x  :=  (xi)f=1  to  indicate  a  point  in  Rn.  Tensors  are  indicated  by  a  further 
underlining:  a  —  ( <Jij)ij=i . 

Suppose  that  the  interior  of  a  compressible  viscoelastic  body  Q  occupies  ft  and  that  its 
surface  coincides  with  5ft.  If  at  a  time  t  this  body  is  subjected  to  a  system  of  body  forces 
/  :=  (/t(iM))£=  i,  for  x  6  ft,  and  surface  tractions  g  :=  (ji(®,t))JLi,  for  x  G  IV,  then  the 
body  Q  will  deform  from  its  equilibrium  configuration.  A  material  particle  originally  at 
the  point  x  will  move  to  the  new  time  dependent  location  x  +  u(x^t)  where  u  :=  (ui)f=1 
denotes  the  displacement  vector.  In  the  linear  theory  these  displacements  define  the 
symmetric  strain  tensor  £  :=  by  the  relationships: 


(dui  duj  \ 
dxj  +  dxi  J  ’ 


(1.4) 


In  addition  to  this  strain  field  there  will  also  be  induced  in  Q  a  stress  field  described  by 
the  symmetric  stress  tensor  er  :=  This  stress  field  rationalizes  the  internal  force 

field  which  is  set  up  within  Q  to  resist  the  external  forces  /  and  g. 

The  stress  field  can  be  related  to  u,  /  and  g  by  Newton’s  second  law  of  motion  (see 
later  in  equation  (4.1))  and  so  it  is  of  interest  to  derive  a  constitutive  relationship 
linking  a  and  u,  or  in  practice,  linking  the  tensors  a  and  £. 
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Constitutive  relationships 


2.1  Hereditary  constitutive  relationships 


In  classical  linear  elasticity  theory  the  constitutive  relationship  between  stress  and  strain 
is  provided  by  Hooke’s  law: 

Oij  =  DijkiEki  or  a  =  D  e, 

where  D  is  a  positive-definite  fourth-order  tensor  of  elastic  coefficients  satisfying  the  sym¬ 
metries 


Dijki  —  Djikh  &ijki  —  Dijiki  and  Dijki  DkUj • 

The  first  two  of  these  are  implied  by  the  symmetry  of  a  and  e  while  the  third  follows  from 
energy  considerations.  However,  in  viscoelasticity  the  third  of  these  only  applies  when  the 
material  is  isotropic,  see  [30,  Equations  (1.10)  and  (2.62)]. 

One  way  of  deriving  a  constitutive  relationship  for  viscoelastic  materials  is  to  assume 
that  a  Boltzmann  superposition  of  stress  increments  can  be  applied  where  these  stress 
increments  are  related  by  Hooke’s  law  to  corresponding  strain  increments.  For  example, 
suppose  that  Q  is  quiescent  for  t  <  0  so  that  e(t)  =  0  for  t  <  0,  and  that  at  t  =  0  the  body 
undergoes  a  strain  c(0).  Then  for  t  >  0  the  resulting  stress  is  assumed  to  be  given  by 

ffoW  =  D(t)e(  0), 

where  a  time  dependence  has  been  introduced  into  the  Hooke’s  tensor  D.  Physically  we 
expect  D  to  be  a  smooth  monotone  decreasing  function  of  t  since  it  is  unrealistic  to  expect 
q  to  grow  over  time  for  the  fixed  strain  e(0).  (Where  would  the  strain  energy  come  from?) 
In  fact  experiments  on  polymers  show  that  D  does  in  fact  decrease  and  this  phenomena 
is  known  as  stress  relaxation. 

Now,  let  At  be  a  small  time  interval  and  set  t{  ~  iAt.  We  approximate  the  strain 
evolution  by  the  step  function 

€{t)  &{ti)  in  [£j,  ti+\)  for  i  0, 1, 2, ... , 
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and  then  each  strain  increment, 


. —  £(^i-fl)  &{ti) , 


induces  a  stress  increment  according  to  Hooke’s  law: 

A Qj{ti)  :=  D(ti  -  tj)Ae(tj)  for  1  <  j  <  i. 

Notice  that  each  of  these  stress  increments  will  also  relax  according  to  the  time  dependence 
of  D.  The  total  stress  at  time  L  is  now  given  by  superposition: 

t 

<r(t»)  :=  <To(*i)  + 

j=  1 

i 

and  by  taking  an  appropriate  limit  we  get  the  hereditary  constitutive  law  as 

a(xyt)  =  D(t)e(u(x,Q))  +  [  D(t- s)e(u'(x,s))ds.  (2.1) 

Jo 

Since  we  are  assuming  that  D(t)  is  smooth  we  can  arrive  at  an  alternate  form  by  partial 
integration, 

<r(aj,  t)  =  D(0)e(u(x ,  £))  —  [  Ds(t  —  s)e(u(x ,  s))  ds ,  (2.2) 

where  the  subscript  s  indicates  partial  differentiation  with  respect  to  the  history  variable 
s.  Either  of  these  may  be  used  as  the  constitutive  relationship,  and  each  demonstrates 
clearly  the  role  of  memory  in  viscoelastic  modelling. 

To  get  a  feel  for  the  form  of  the  time  dependence  of  the  stress  relaxation  tensor  D  we 
now  describe  a  perhaps  more  intuitive  method  for  deriving  these  constitutive  relationships. 


2.2  Spring  and  dashpot  models 

We  start  with  the  physical  observation  that  viscoelastic  materials  display  the  characteris¬ 
tics  of  both  elastic  solids  and  viscous  fluids.  The  kinetics  of  these  type  of  substances  are 
modelled  respectively  by  the  spring  and  the  dashpot.  In  these  models  the  stress  carried 
by  the  spring  is  proportional  to  the  strain  in  the  spring  and  is  given  by  Hooke’s  law: 
a  =  Ee .  The  stress  carried  in  the  dashpot  is  proportional  to  the  strain  rate  and  is  given 
by  Newton’s  law  of  viscosity:  a  =  rje*. 

One  then  models  a  viscoelastic  material  by  considering  a  notional  system  of  springs 
and  dashpots  with  independent  stiffness  and  viscosity  parameters.  There  are  essentially 
two  ways  to  connect  a  spring  to  a  dashpot:  in  series  and  in  parallel.  These  are  the 
building  blocks  and  are  named  the  “Maxwell”  and  “Voigt”  models. 


2.2.  SPRING  AND  DASHPOT  MODELS 
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Figure  2.1:  A  HOOKEAN  (linear)  SPRING:  a  =  Ee\  E  IS  THE  SPRING  STIFFNESS 


E 

AAAA/V 


e,  a  =  Ee 


de 

Figure  2.2:  A  Newtonian  (linear)  dashpot:  a  =  rj— ;  77  is  the  viscosity 


V 


e,  o  =  r) 


de 

dt 


The  Maxwell  model 


Figure  2.3:  The  Maxwell  model 


E 


I  £5, 0-5 


AAAAAr 


*1  |  ££>,  CT£) 

3—^ 


The  Maxwell  model  is  a  series  connection  of  a  spring  and  dashpot.  In  this  model  es  and 
as  denote  the  strain  and  stress  in  the  spring  alone,  and  en,  &d  denote  those  in  the  dashpot 
alone.  The  total  stress  is  given  by  a  —  as  —  cfd  and  the  total  strain  by  e  =  Es  +  £d- 
Differentiating  and  using  Hooke’s  and  Newton’s  laws  yield 

=  +  ^  +  (2  3) 

dt  E  dt  rj  dt  t  dt'  v  '  ’ 

where  r  :=  t]/E  is  the  so-called  relaxation  time.  Using  <r(0)  =  Ee(0)  this  ODE  is  easily 
solved  to  give 

a(t)  =  Ee-VTs{0)  +  E  f*  e-^'T£'(s)  ds, 

Jo 

and  this  is  essentially  (2.1)  with  the  scalar  analogue  of  D  given  by  D(t)  =  Ee~t/T. 


The  Voigt  model 


Connecting  the  spring  and  dashpot  in  parallel  yields  the  Voigt  model.  This  time  ss  = 
ed  =  e  and  equilibrium  demands  that  a  —  as  +  od,  hence 


de 

"dt+Ee 


de  e 


a 

V* 
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Figure  2.4:  THE  VOIGT  MODEL 


This  gives  the  constitutive  law  in  hereditary  form  as 


e(t)  =  e-‘/Te(0)  +  -  [' e-^s)/Ta(s)ds. 
V  Jo 


The  Maxwell  solid 


In  his  internal  variable  formulation  A.R.  Johnson,  in  for  example  [24],  uses  these  basic 
building  blocks  in  the  Maxwell  solid.  Here  E0  and  Ey  are  spring  stiffnesses  and  a*,  e*  are 
internal  stress  and  strain  variables.  This  time  a *  —  Eye* ,  bd  =  e  —  e*  and  as  =  EqBs- 
Also  a*  =  <jd  and  this  gives 


Eie*=r1jt(e-e*) 


de*  e*  _  de 
dt  +  t  dt ’ 


where  now  r  :=  tj/E\.  Solving  this  we  get 

e*(t )  =  e-t/Te(0)  +  /  e~^~s^Te'{s)  ds  =  e{t)  —  —  f  e~^~s^Te(s)ds.  (2.4) 

Jo  T  Jo 

where  e(0)  could  be  any  constant  but  arises  here  from  the  initial  condition  e*(0)  =  e(0). 
Now,  defining  the  stress  relaxation  function 


D(t)  :=  E0  +  Eie-^T 


as  the  scalar  analogue  to  the  tensor  D(t)  in  (2.1)  and  (2.2),  and  using  this  in  (2.4)  along 
with  the  relation 


a  -  as  +  a*  -  E0e  +  Eye*  (since  es  =  e), 


gives 

a(t) 


E0e(t)  +  Eie  t/Te(0)  +  J  Eye  (t  s^Te'(s)ds, 
D(0)e(t)  -  [  Ds(t  -  s)e(s)ds. 
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Figure  2.5:  The  Maxwell  solid 


£,  a 


This  is  the  scalar  analogue  of  equation  (2.2)  and  suggests  that  we  model  D  with  the 
Dirichlet-Prony  series, 

D(t)=cp(t)D(0)  (2.5) 

where  (p(t)  is  a  generic  stress  relaxation  function  given  by 
N 

cp(t)  —  (po  +  (2*6) 

t=l 

Here  the  (possibly  x  dependent)  coefficients  {<Pi}iLo  316  non- negative  and  normalized  so 
that  <p(0)  =  1,  and  the  (possibly  x  dependent)  {<*{}■) Lx  are  non-negative.  More  generally 
one  could  of  course  write  (with  summation  not  implied), 

-W*  jkl 

Dijkl{t)  :==  {Dijkl) 0  +  (Dijkijm  exP("”(CKijfc/)m^)* 
m— 1 

The  Dirichlet-Prony  series  is  an  extremely  convenient  form  to  take  for  large  scale  compu¬ 
tational  approximations  to  problems  (1.1),  (1.2)  and  (1.3)  since  if 

m  ~ 

then  one  can  exploit  the  simple  recurrence 
■4>{t  +  A:)  =  e-QV(i) 

to  update  the  history  term  arising  from  a  discretization  of  the  Volterra  integral.  For  general 
Volterra  problems  one  must  usually  store  the  entire  solution  history  as  the  computation 
advances  through  the  time  levels  and  moreover,  at  each  time  level  this  history  needs  to  be 
summed  to  approximate  the  integral.  For  such  methods  the  number  of  operations  required 
at  time  level  N  is  of  the  order  0(N2).  The  Dirichlet-Prony  series  provides  a  very  useful 
short  cut  around  this  UN2  problem” .  (In  certain  special  cases  one  can  also  overcome  this 
difficulty  using  other  means,  see  for  example  [22,  19]). 
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2.3  The  Reversed  Maxwell  Solid 

One  can  also  arrive  at  a  Maxwell  solid  by  switching  the  order  in  which  the  spring  and 
dashpot  appear  in  the  viscous  arm  of  the  network.  We  call  this  the  Reversed  Maxwell 
Solid. 


Figure  2.6:  The  Reversed  Maxwell  Solid 
V  I  <?■*,£* 


-WAV 


E0 

-WAV 


e,  a 


Balancing  stresses  such  that  a *  =  a1  and  using  el  =  e  -  e*  now  gives  the  differential 
equation  for  the  new  internal  variable  e*  as, 

^ -  -j-  —  —  — ,  where  r  :=  rj/Ex. 
at  r  r 

Note  that  £*  is  quite  different  from  e*  as  introduced  earlier  for  the  Maxwell  model,  although 
there  is  a  simple  connection  which  we  demonstrate  below.  This  time  we  have, 

e*{t)  =  -  f  e~^~s^re{s)ds, 

T  Jo 

where  we  used  £*(0)  =  0  (otherwise  any  term  of  the  form  Ae~t/T  may  be  added  on.). 
Since  es  =  e,  the  total  stress  a  is  given  by, 

a(t)  =  as  +  cr1  =  EqSs  +  E\El , 

—  Eq£  +  Ei(e  —  £*), 

=  (So  +  EMt)  -  ^  T  e-(*-)/Te(«)  ds, 

T  JO 

=  D(0)e(t)  -  [*  Ds(t-s)e(s)ds, 

Jo 

where  again  we  set, 

D(t)  :=  S0  +  Sie-‘/T. 

This  is  exactly  as  before  for  the  Maxwell  solid,  and  so  again  represents  a  scalar  analogue 
of  (2.2).  The  question  is  whether  or  not  this  reversed  model  is  of  any  use,  and  in  this 
context  we  note  the  possibility  of  strong  stability  estimates  for  solution  of  this  ODE. 


2.4.  THE  GENERALIZED  MAXWELL  SOLID:  INTERNAL  VARIABLES 
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Note  that  e*  and  e*  are  very  simply  related  by, 


e(t)  —  e*(t)  +  £*(f). 


This  is  obvious  since  the  total  strain  in  the  viscous  arm  of  the  networks  is  the  sum  of  the 
strains  in  each  component.  This  relationship  is  easily  proven  by  considering  the  integral 
representations  of  e*{t)  and  £*(f)  given  earlier. 

Let  us  now  briefly  try  to  make  a  connection  with  the  “ODE  formulation”  as  described 
by  Janovsky  et  al.  in  [22].  Using  e  =  e*  +  e*  in  the  differential  equation  for  e,  we  easily 
obtain  a  differential  relationship  between  e*  and  £*  as, 

cfe*  _  e* 
dt  t 

The  internal  variable  £*  plays  essentially  the  same  role  as  (the  scalar  analogue  of)  w  in 
[22,  Equation  4.1]. 

Finally  in  this  section  we  note  that  we  can  also  arrive  at  ODE’s  for  the  internal  stress 
variables.  For  the  Maxwell  solid  we  have, 


-jr  +  —  =  EX  — 
dt  r  dt 

and  for  the  Reversed  Maxwell  solid  we  get: 

TTT  +  E\£*  =  E\e, 

at 

do*  de*  de  _  do*  o* 

=>  -dr+El~dF-El'dt-~dr+ v 

Clearly  o *  =  o*  and  so  (again)  de*/dt  =  e*/r. 


2.4  The  generalized  Maxwell  solid:  internal  variables 

We  now  return  to  the  Maxwell  solid  and  generalize  the  conceptual  spring  and  dashpot 
model  in  order  to  motivate  the  choice  of  the  Dirichlet-Prony  series  for  the  relaxation 
function  as  given  in  (2.6).  To  begin  with  we  assume  again  a  state  of  uniaxial  stress  and 
strain. 

The  generalized  Maxwell  solid,  shown  in  Figure  2.7,  consists  of  a  Hookean  spring 
connected  in  parallel  to  a  sequence  of  N  spring-dashpot  components.  In  this  model, 

£o  =  e,  cro  =  and  o*  =  E^*. 

Balancing  the  stresses  carried  by  each  of  the  spring-dashpot  pairs  we  get  for  each  i  £ 
{1, . . .  ,N}  that, 

de*  e*  _  de 
dt  +  T{  dt ’ 

=>  eUt)  =  e-t/T<£(0)  +  e"(‘-a)/V(s)  ds  =  e(t)  -  -  [*  e-(*-*)Me(a)  da, 

Jo  Ti  Jo 
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Figure  2.7:  The  Generalized  Maxwell  solid. 


where  we  used  s*(0)  =  e(0)  and  this  time  we  have  set  :=  Ei/r]i.  The  total  stress  carried 
by  the  assemblage  is  therefore  given  by: 

a(t)  =  <ro(<)  +  <7i (t)  H - t-0iv(t), 

=  Eoe(t)  +  E\e\{t)  H - f  Eff£*N(t), 

=  E0e(0)  +  E0(e(t)  -  e(0)) 

N  ,  rt  x 

+  E  (^e-t/T<e(0)  +  Jo  Eie-^-^e'is)  dsj  , 

=  E{t)e{ 0)  +  ^  E(t  -  s)e'(s)  ds,  (2.7) 

Jo 

where 

N 

E{t)  -^Eo  +  ^Eie-^. 

i=l 

The  constitutive  relationship  (2.7)  is  the  scalar  analogue  of  (2.1)  with  the  analogue  of 
D(t  —  s)  given  by  E(t  —  5),  which  itself  is  an  example  of  the  Dirichlet-Prony  series  given 
in  (2.6).  Note  that  if  we  set  Eq  :=  0  then  this  generalized  Maxwell  solid  actually  models 
a  fluid  since  YimE(t)  =  0. 


2.5.  THE  GENERALIZED  REVERSED  MAXWELL  SOLID 
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2.5  The  generalized  Reversed  Maxwell  Solid 

By  switching  the  order  of  the  springs  and  dashpots  in  each  of  the  viscous  arms  of  the 
generalized  Maxwell  solid  network  we  arrive  at  the  generalized  Reversed  Maxwell  Solid. 

Figure  2.8:  The  Generalized  Reversed  Maxwell  Solid. 


This  time  we  balance  the  stresses  in  each  spring-dashpot  pair  such  that  a\  =  E{{e— el) 
and  obtain, 


ae;  e;  _  e 

dt  Ti  Ti 


for  each 


and  where  t*  :=  rji/Ei.  Again,  with  £^(0)  =0  we  can  solve  for  the  internal  strains  to  get, 
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as  expected,  and  we  obtain  the  total  stress  from: 

a(t)  =  Co  +  E\(e  —  ei)  4 - +  Eff(e  —  e^), 

=  (E0  +  •  •  •  +  EN)e(t)  -  [l  +  •  •  •  +  e(s)  ds, 

Jo  V  n  tat  / 

=  £?(0)e(^)  -  [  Es(t  -  s)e(s)  ds , 
where,  again,  we  set 
N 

E(t)  :=E0  +  J2Eie-t/Ti. 

i— 1 


This  constitutive  relationship  between  a  and  e  is  a  scalar  analogue  of  (2.2),  which  itself 
is  equivalent  to  (2.1).  Hence  the  constitutive  relationship  generated  by  the  generalized 
Reversed  Maxwell  Solid  is  equivalent  to  that  generated  by  the  Maxwell  Solid. 

Once  again  it  is  trivial  to  note  the  connection, 


e(£)  —  (t)  +  el(t ), 


for  each  i,  and  it  follows  that, 


d4  =  e[ 

dt  Ti  * 


Introducing  the  internal  stresses  we  also  get, 

d4 

dt 


rji— f  +  Eie\  —  EiS , 


daj  del  _Fde  _  dai  ■  <?* 
dt  1  dt  ldt  dt  Ti’ 


<r*  —  cr*  =  constant  =  0, 


since  the  stress  in  each  arm  is  independent  of  the  order  in  which  the  components  are 
arranged. 


2.6  Other  constitutive  relationships 

The  Dirichlet-Prony  series  is  not  however  the  only  form  used  to  model  the  stress  relaxation 
functions,  for  example  the  authors  of  [1]  use  the  stretched  relaxation  function 

<p(t)  =  (fi0  exp (~(at)p)  for  p  £  (0, 1].  (2.8) 

Obviously  no  simple  recurrence  exists  for  this  form.  Another  popular  choice  for  <p  is  the 
power  law  where 


<p(t)  =  (pot  v 


forpe  (0,1), 


(2.9) 


2.6.  OTHER  CONSTITUTIVE  RELATIONSHIPS 
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although  from  either  of  (2.1)  or  (2.2)  this  implies  that  either  e(0)  is  zero  irrespective  of 
the  magnitude  of  the  load,  or  g-(O)  is  infinite.  Neither  of  these  are  physically  realistic  and 
so  we  would  prefer  to  modify  this  law  to 

<p(t)  =  tpo(t  +  <pi)~p  for  pe  (0,1),  (2.10) 


where  (pi  >  0  in  order  to  remove  the  non-physical  behaviour.  Nonetheless,  it  is  instructive 
to  see  how  one  might  “derive”  the  power  law,  and  for  this  we  borrow  heavily  from  Chern’s 
thesis  [4]  which  exploits  the  fractional  calculus. 

The  formulation  is  based  on  the  observed  fact  that  viscoelastic  materials  behave 
in  a  way  intermediate  to  that  of  solids  and  fluids.  Interpreting  this  literally  yields  a 
constitutive  law  that  contains  fractional  derivatives.  Unfortunately  we  are  unable  here  to 
give  this  interpretation  the  depth  it  deserves  and  instead  try  only  to  illustrate  the  main 
point.  Recall  that  the  stress  in  a  solid  is  proportional  to  the  strain  while  the  stress  in  a 
fluid  is  proportional  to  the  strain  rate.  Accepting  the  intermediate  nature  of  viscoelastic 
materials  the  idea  is  to  define  the  viscoelastic  constitutive  law  as: 

a{t)  =  D{0h(t)  +  £{1)0?e(f),  (2.11) 

for  constant  fourth  order  tensors  and  D^\  and  where  a  €  [0,1).  The  fractional 
derivative  operator  may  be  defined  as: 

dt£(t)  ~  (fp- Z  a)  JQ  (*  “  s)'“e(s) ds)  .  for  a  €  [0, 1).  (2.12) 


(Note  that  a  can  be  irrational,  even  though  the  word  “fractional”  is  always  used.)  By  firstly 
integrating  by  parts  in  (2.12)  and  then  taking  the  differentiation  through,  Chern  arrives  at 
a  constitutive  law  which  is  suitable  for  use  within  the  standard  finite  element  framework. 
Two  solution  schemes  are  considered:  a  solution  in  the  Laplace  transform  domain  and  a 
direct  time  domain  solution.  However,  in  this  case  there  is  no  efficient  history  storage  and 
so  the  operation  counts  and  computer  memory  requirement  grow  without  bound  as  the 
time  step  is  diminished. 

The  “justification”  for  the  power  law  is  as  follows.  Carrying  out  this  integration- 
differentiation  process  gives 


d?e(t)  = 


r(i 


-a)£{0)  +  F(T 


hr)H{t-s)~a£,{s)ds ’ 


(2.13) 


and  using  this  in  the  scalar  analogue  of  equation  (2.11)  we  now  arrive  at  the  constitutive 
law: 


f  ®  JpJ-t  H 

a(t)  =  E0e(t)  +  r/1  _  ^  6(0)  +  r(1_a)  fQ(t-  s)'V(s)  ds. 


T(l-a) 


(2.14) 


This  seems  to  combine  (2.1)  and  (2.2)  when  (p(t)  is  given  by  the  power  law,  (2.9). 

We  now  have  several  candidates  for  the  constitutive  law  and  these  may  be  used  to 
generate  a  variety  of  differential  equation  problems.  Later  in  Chapter  4  we  do  just  this 
and  demonstrate  how  concrete  forms  of  the  abstract  problems  (1.1),  (1.2)  and  (1.3),  as 
well  as  some  non-standard  variants,  can  be  derived  to  model  viscoelastic  behaviour. 


Chapter  3 


Multidimensional  internal 
variables 

3.1  Overview 


So  much  for  uniaxial  states  of  stress  and  strain.  In  fact  it  can  be  shown  that  for  each 
relaxation  mode  (i.e.  each  spring-dashpot  pair)  in  higher  dimensions,  there  is  an  ODE 
governing  the  evolution  of  each  of  the  internal  strain  tensor  components.  Thus  we  have, 

J  *rrn  rfcr.  . 

- -I - lJ-  =  — 11  (for  fixed  x), 

dt  7 ~<ffi  dt 

where  the  details  are  given  below.  The  significance  of  these  internal  variable  formula¬ 
tions  for  the  viscoelastic  constitutive  behaviour  lies  in  the  fact  that  it  is  possible  to  solve 
some  kinds  of  viscoelasticity  problems,  when  the  relaxation  functions  are  in  the  form  of  a 
Dirichlet-Prony  series  (2.6),  using  only  a  linear  elasticity  solver  and  an  ODE  solver.  This 
obviates  the  need  to  create  special  software  for  quasistatic  viscoelasticity  problems.  For 
more  on  this  we  refer  again  to  [24]  and  also  to  [43],  but  we  return  to  this  topic  in  a  later 
chapter. 


3.2  Multidimensional  constitutive  relationships 


In  [54]  for  example  Shaw  and  Whiteman  use  a  hereditary  multidimensional  constitutive 
relationship  of  the  form, 

Oij{t)  =  Dijkl(0)eki(t)  -  dDl^ — —eki{s)  ds 

at  each  fixed  point  a;  in  a  viscoelastic  body  (see  (2.2)  given  earlier).  Here  D(t)  = 
k  l-\  (f°r  Problems  posed  in  Rn)  is  a  fourth-order  stress  relaxation  tensor 


18 


3.2.  MULTIDIMENSIONAL  CONSTITUTIVE  RELATIONSHIPS 


19 


given,  for  example,  by  the  Dirichlet-Prony  series, 

No 

Dijkl{t)  -=  {Dijkl) 0  +  ^  ( Dijkl)me  ^  m- 

m— 1 

However,  Johnson  and  Tessler  in  [24,  equation  20.9]  write  this  in  a  different  form  and  it 
is  useful  here  to  explore  the  connection. 

Assuming  that  t  =  0  is  a  reference  time  such  that  e(t)  =  0  for  all  t  <  0,  then  the 
constitutive  equation  is, 

a„(t)  =  Cijkieul.t)  +  J‘_^ 

=  {Cm+  -CW 0))et,(i)  -  f‘  °  ~  S)eu(s)  is. 

By  comparing  these  two  forms  we  see  that 


Dijki{0) 

II 

£ 

£ 

+ 

*  Cijki  (0) , 

and  also, 

Dijklit) 

=  *C'ijkl(t), 

Dijki{t) 

—  *Cijki(t)  +  constant;^, 

Dijkl  (t) 

=  *cijki(t)  + Cijki  Vt. 

In  fact  Johnson  and  Tessler  in  [24,  equation  20.7]  write, 

Nd 

:=  £  C:-Tie-t/Tm  for  t  >  0, 
m=  1 

where  C*-\i,  C*^ , . . . ,  C*j^D  are  temporally  constant  fourth-order  tensors.  (In  fact  this  is 
a  slight  generalization  of  Johnson  and  Tessler’s  expression,  but  no  matter.) 

Comparing  terms  once  again  with  our  Dirichlet-Prony  series  representation  of  D  we 
see  that  simultaneously  we  must  have, 

Nd 

Dm{t)  =  C'iM+ZCft*-*'", 

m= 1 
ND 

Dijkl(t)  =  (Ajw)o+  E(Ajfci)me“t/%, 

771=1 

and  so  we  infer  that, 

(Dijki)  o  =  Cijki  and 


for  each  m. 
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CHAPTER  3.  MULTIDIMENSIONAL  INTERNAL  VARIABLES 


3.3  Multidimensional  internal  variables 

It  is  of  interest  to  derive  “ODE”  representations  of  these  hereditary  constitutive  relation¬ 
ships  in  terms  of  internal  strain  variables.  To  begin  we  set, 

:=  e~t/rmeki{ 0)  +  [  e~^~s)/Tm  de^S)-  ds  =  eki(t )  -  —  f  e~^~s),Tmeki{s)  ds, 
Jo  os  Tm  Jo 

and  then  by  differentiating  we  obtain, 

d.:±m  =  _  J_  [e_t/rmew(0)  +  t  e~(t~s)/Tm^^-ds]  . 

dt  dt  rm  Jo  ds  J 

Substituting  for  the  integral  now  yields  the  family  of  differential  equations  governing  the 
evolution  of  the  internal  strain  variables: 

a  'em  .  *e£K<)  deki(t) 

dt  rm  dt 

For  the  stresses  we  note  that, 

Nd 

Oi j(t)  =  Cijki£ki(t)  +  ^2  C*jki  *£™i(t)- 
m=I 

We  can  also  derive  ODE’s  for  the  internal  stresses  defined  by, 

:=  C*$ 

which  gives, 

nd 

CTijit)  =  Cijkl£kl(t)  +  J2  *<$(*)• 

171=1 


In  this  case  the  evolution  equations  are, 

a  VW(t)  _  *>mdekl(t) 

dt  +  rm  ijkl  dt  ' 

3.4  Reversed  multidimensional  internal  variables 

From  our  consideration  of  the  Maxwell  models  we  define  the  “reversed”  internal  strain 
tensors  *£m(t)  via, 


•eSW  :=  eu(t)  -  *ej5W, 

and  of  course  it  follows  that, 

.e3(<)  =  —  f  e-(t-a)/Tmeki(s)ds. 

Tm  Jo 
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Then, 

dt  Tffi 

similar  to  before,  and  also, 

dt  Tm  Tm 

These  are  clear  analogues  of  the  results  for  the  scalar  case  given  earlier. 


Chapter  4 


Model  problems  in  viscoelasticity: 
Volterra  formulations 


4.1  Viscodynamics 


To  obtain  the  governing  equations  for  the  dynamic  response  of  a  viscoelastic  body  one  uses 
Newton’s  second  law  to  relate  the  stress  field  a  and  the  forces  f  and  g  to  the  acceleration, 
or  inertia,  of  the  body  Q.  This  process  is  familiar  from  linear  elasticity  theory  and  gives, 
with  boundary  and  initial  data,  the  following.  For  i  =  1, . . . ,  n: 


'/  —  rr- 
'i  °iJ,3 

=  fi 

in  fixl, 

Ui 

=  0 

in  rD  x  X, 

(JijTlj 

=  9i 

in  Fn  x  X, 

Ui{x ,  0) 

—  Uio 

in  fi, 

ui(xi  0) 

-  Uil 

in  fi. 

> 

Here:  repeated  indices  imply  summation;  X  (0 ,T)  is  a  time  interval;  g  is  the  mass- 

density  of  Q\  and,  ft  (ni)”=1  is  the  unit  outward  directed  normal  to 

Using  (2.2)  to  substitute  for  the  stress  one  arrives  at  the  Partial  Differential  Volterra 
(pdv)  problem:  find  u  such  that 

Qu'lit)  -  (Dijkl(0)ekl(umtj  =  Mt)  -  j*  (dDljk£~  s)ekl(u(s)))  ds, 


in  fi  x  X  with  the  indicated  initial-boundary  data.  With  an  appropriate  definition  of  A 
and  2?(t,s)  this  is  clearly  a  realization  of  the  abstract  problem  (1.3).  Note  that  it  is  “safe” 
to  use  the  Dirichlet-Prony  series  (2.6)  or  modified  power  law  (2.10)  in  this  problem,  but 
we  may  not  use  the  power  law  (2.9)  directly  because  we  cannot  then  interpret  £>(0). 

In  terms  of  existence,  uniqueness  and  stability  of  solutions  this  problem  has  been 
studied  in  [10,  11,  34].  Numerical  schemes  are  given  in  [13,  60,  39,  42]. 
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One  could  also  use  the  fractional  calculus  model  to  substitute  for  a  in  Newton’s 
second  law.  This  will  yield  a  PDV  equation  of  the  form 

£)(])  o 

Qu"{t)  -  {D\jheki(u(t)))  j  =  fi(t)  +  r J0  ^  ~  s)_“£fc«(«(s))  ds. 

On  the  other  hand  one  could  use  (2.1)  and  then  arrive  at 

gu"(t)  —  (Dijki{t)eki{u(0)))  j  =  fi{t)  +  f  ( Dijki(t  -  s)eki{u' (s)))  j  ds. 

JO 

Note  that  u  does  not  occur  as  a  natural  “unknown”  in  this  problem  and  so  it  is  possible 
to  replace  u  with  u'  and  arrive  at  the  alternative  problem:  find  u  such  that 

rt 

gu[(t)  +  /  (Dijki{t  —  s)£w(w(s)))  j  ds  =  fi(t)  —  (Dijki(t)eki(uo))  j, 

J  0 

which  makes  sense  if  tio  is  smooth  enough.  The  initial  datum  for  this  problem  is  now 
•u(O)  =  tti.  Properties  of  the  solution  of  these  type  of  problems  are  studied  in  [11,  34]  and 
numerical  analysis  is  given  in  [35,  32]. 

However,  one  must  resist  the  temptation  to  interpret  this  as  a  parabolic  problem  for, 
in  general,  it  is  not.  To  see  this  use  the  power  law  (2.9)  with  (2.5)  to  obtain  (with  g  =  1 
and  D  not  *  dependent  for  simplicity): 

u'i(t)  +  Dijki  f  (t-s)~p(eki(u(s))jds  =  fi{t),  (4.2) 

J  0 

where  J  now  incorporates  the  additional  term  in  uo -  In  the  case  p  =  |  we  find  that  the 
operator  I  defined  by, 

1  f*  _  i 

Iw(t)  :=—=  I  (t-s)  2  w(s)  ds 
yTT  Jo 

has  the  property, 

I2w{t)  =  I(Iw){t)  =  [tw(s)ds, 

Jo 

and  so  may  be  regarded  as  the  square  root  of  the  definite  integral  operator.  Applying 
1  1 
d(2  to  both  sides  of  (4.2)  in  the  case  p  =  5  we  arrive  at 

( 2  ««(<)  +  VxDijki(m(u(t)),j  =  dim. 

This  equation  is  half  way  between  being  parabolic  and  hyperbolic.  Similar  manipulations 
are  also  possible  in  the  case  p  ^  with  the  final  time  derivative  being  of  order  between 
1  and  2.  Numerical  methods  for  fractional  order  differential  equations  are  studied  in 
[41,  31,  12]. 

For  more  detail  on  these  type  of  problems  see  [37],  as  well  as  the  other  papers  in  that 
collection. 
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4.2  Viscostatics 

Recall  that  the  classical  linear  elasticity  equations  are  “derived”  from  Newton’s  second 
law  (4.1)  by  dropping  the  inertia  term  guff.  This  corresponds  to  modelling  the  problem 
at  times  long  after  the  load  has  been  applied  when  the  transient  response  has  died  out, 
and  results  in  a  very  well-known  elliptic  problem.  A  similar  approach  can  be  adopted 
for  viscoelastic  response  although  this  time  it  is  a  true  approximation  since  the  resulting 
problem  is  not  time  independent  due  to  the  persistence  of  the  Volterra  term.  It  seems 
that  this  approximation  can  be  useful  when:  the  inertia  term  is  negligible,  which  may 
occur  when  the  load  is  smoothly  and  slowly  applied,  or  even  when  the  applied  load  is  such 
that  the  dynamics  of  the  response  is  dwarfed  by  the  forced  response;  or  when  it  is  the 
long-time  creep  response  that  is  of  interest.  Since  the  time  dependence  persists  we  refer 
to  the  resulting  problems  as  modelling  quasistatic  viscoelastic  response. 

The  governing  equations  for  these  type  of  problem  are  obtained  from  (4.1)  by  setting 
Qun(t)  :=  0,  setting  J  :=  [0,T]  and  discarding  the  initial  data.  Thus,  for  i  =  1, . . .  ,n  we 
have, 


&ij,j  —  fi  in  f2  x  Cf , 

U{  =  0  inFjoxJ', 

crijrij  —  Qi  in  TN  x  {J , 

which  are  turned  into  differential  equation  problems  for  u  by  substituting  for  the  stress 
using  either  of  (2.1)  or  (2.2).  These  give  respectively  the  PDV  problems:  find  u  such  that 
for  each  i  €  {1, .  - . ,  n}, 

~  —  s)eki{u  (5)))  j  ds  =  fi(t)  +  (Dijki(t)eki(u( 0))) 

and 

-(Dijki(0)eki(u(t)))tj  =  fi(t )  -  (dD'3Ufo — —*«(«(*)))  ds. 

The  first  of  these  is  essentially  a  Volterra  first-kind  equation  for  u\  while  the  second  is  a 
second-kind  equation  for  u.  In  both  cases  one  obtains  u(0)  by  solving  a  linear  elasticity 
problem. 

Numerical  schemes  and  a  priori  error  estimates  were  first  provided  for  both  of  these 
problems  in  [45].  Later  and  for  the  second-kind  problem  only,  the  estimates  were  improved 
(in  terms  of  the  size  of  the  error  constant)  in  [44].  These  latter  results  depend  on  by-passing 
Gronwall’s  inequality  and  using  more  sensitive  comparison  results  to  obtain  sharp  data- 
stability  estimates.  These  estimates  have  now  been  generalized  in  [57].  Also  for  the  second- 
kind  problem,  a  posteriori  error  estimates  for  a  space-time  finite  element  discretization  of 
a  model  problem  have  been  given  in  [48]  and  [52].  These  results  are  based  on  the  error 
estimates  in  [47]  and  are  currently  being  generalized  to  the  multidimensional  problem 
described  above  in  [51],  [54],  [55]  and  [56]. 


4.3.  NON-FICKIAN  DIFFUSION 
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4.3  Non-Fickian  diffusion 

This  section  is  not  directly  relevant  to  the  Seed  Project,  but  it  is  another  illustration 
of  where  viscoelastic  effects  play  an  important  role,  and  therefore  need  to  be  modelled 
accurately.  We  include  this  material  for  interest  only. 

In  classical  diffusion  theory  the  gradient  of  the  concentration  u  of  an  active  agent  (the 
penetrant)  diffusing  through  a  carrier  medium  is  related  to  the  mass  flux  by  Fick’s  law: 
J  =  —A Vu,  where  A  is  the  diffusivity  of  the  carrier  substance.  Conservation  of  mass  then 
demands  that  u '  =  —  V  •  J  which  yields  the  familiar  heat  equation, 

v!{t)  =  V  •  AVu. 

If  we  define  M(t)  as  the  total  mass  of  penetrant  absorbed  by  the  carrier  per  unit  area 
at  time  t  then  it  is  well  known  (from  similarity  solutions)  that  M(t)  ~  t*  for  Fickian 
diffusion. 

Diffusion  in  rubbery  polymers,  those  well  above  their  glass  transition  temperature 
(GTT),  is  according  to  Durning  in  [14]  adequately  described  by  Fick’s  law,  but  the  sit¬ 
uation  is  much  more  complicated  for  glassy  polymers,  those  near  but  above  their  GTT. 
As  the  penetrant  moves  through  the  polymer  it  can  force  a  phase  change  and  so  leave 
behind  it  the  polymer  carrier  in  its  rubbery  state.  The  stiffness  and  relaxation  properties 
of  the  polymer  change  abruptly  by  orders  of  magnitude  across  this  phase  change  (see  for 
example  [18]),  and  as  a  result  a  differential  stress  is  set  up  across  the  penetrant  boundary. 
Moreover,  because  the  carrier  is  viscoelastic  this  stress  is  described  by  a  hereditary  consti¬ 
tutive  law  and  this  behaviour  provides  a  mechanism  for  the  observed  non-Fickian  effects. 
Workers  in  the  field  make  the  following  very  rough  classification. 

Case  I  diffusion:  standard  Fickian  diffusion  where  M(t)  ~  1 2,  applies  to  polymers  in 
the  rubbery  state  high  above  the  GTT. 

Case  II  diffusion:  non-Fickian  diffusion,  M(t)  ~  ta  where  \  <a<  1,  applies  to  glassy 
polymers  near  to  but  above  the  GTT. 


There  is  also  a  “Super  Case  II”  category  corresponding  to  a  >  1,  see  [6].  For  Case  II 
sharp  fronts  (rather  like  shocks)  may  appear  as  the  penetrant  diffuses  through  the  carrier. 
This  front  moves  initially  at  a  constant  speed  and  then  slows  down,  [8],  and  this  explains 
why  M(t)  is  almost  linear,  and  thus  Mf(t) — the  rate  of  absorbtion — is  almost  constant. 
By  contrast  M'(t)  for  Case  I  is,  in  the  words  of  Cox  in  [8],  “delt a- funct ion-like” ,  and  this 
property  of  glassy  polymers  has  an  interesting  application  in  the  area  of  controlled  drug 
delivery  products.  Cox  gives  a  nice  example. 

An  active  agent  (the  drug)  is  embedded  into  a  polymer  through  which  it  cannot 
diffuse.  This  may  for  example  be  a  tablet  which  is  to  be  swallowed.  When  the  carrier  is 
invaded  by  a  solvent,  such  as  digestive  fluid,  the  drug  can  then  diffuse  out  of  the  polymer 
through  the  solvent  in  a  non-Fickian  way.  Since  Mf(t)  is  almost  constant,  this  allows  a 
controlled,  constant-rate  delivery  of  the  drug  to  the  body  for  several  hours. 
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The  polymer  doesn’t  have  to  be  a  tablet.  In  fact,  according  to  Cohen  and  White  in 
[6]  (who  also  describe  other  applications  of  non-Fickian  diffusion),  such  “smart”  pharma¬ 
ceutical  products  can  be  designed  to  be  “swallowed,  smelled,  surgically  implanted,  rubbed 
on,  taped  on,  strapped  on” ,  and  can  in  effect  be  applied  to  any  part  of  the  body.  There  is 
an  extensive  literature  on  this  science  and  in  addition  to  those  already  cited  we  refer  also 
to  [20,  7,  15]. 

To  get  a  flavour  of  the  mathematical  modelling  that  these  researchers  employ  we 
borrow  from  [5]  and  consider  the  modelling  of  one-space  dimensional  diffusion  through  a 
glassy  polymer.  Our  development  yields  a  linear  model,  but  it  is  unlikely  that  this  will 
reproduce  the  sharp  fronts  characteristic  of  polymer  diffusion.  The  references  cited  deal 
with  realistic  nonlinear  models. 

To  account  for  the  differential  stress  set  up  at  the  penetrant  front  Fick’s  law  is  modified 
to  include  a  stress  dependence  in  the  following  way: 


J  —  -(A  ux  +  Kcrx). 

Here  u  is  the  concentration,  A  the  usual  (Fickian)  diffusion  constant,  and  k  is  a  propor¬ 
tionality  constant.  Conservation  of  mass  again  demands  that  uf  =  —Jx  and  this  gives 

u  =  Ait^x  +  ftaxx. 


The  stress  is  viscoelastic  and  the  usual  approach  is  to  adopt  the  Maxwell  model,  given 
earlier  in  (2.3),  with  the  assumption  that  u  depends  linearly  on  strain  rate  e*  (this  is  in 
order  to  get  true  Case  II  behaviour — see  [9]).  Thus, 


da  a 
—  +  -  -  /iU, 


where  fi  is  a  proportionality  constant.  In  the  nonlinear  theory  the  dependence  of  r  on  u 
is  crucial,  but  here  we  shall  assume  that  r  is  constant.  Integrating  we  get, 


<r(t)  =  Me-t/Tu(0)  +  n  [*  e~^~s^Tu(s)  ds. 

Jo 

Eliminating  the  stress  from  the  transport  equation  and  using  mass  conservation  gives  the 
single  differential- Volterra  equation, 

/*£ 

uf(t)  =  A  uxx  +  njie~^ruxx{  0)  +K/J,  eT^~s^Tuxx(s)  ds. 

Jo 


Assuming  for  simplicity  that  u(0)  =  0  we  can  generalize  this  to  a  multidimensional  model 
and  obtain  the  PDV  equation, 

u'(t)  =  V  •  A Vu  +  V  •  (/sV  J*  ne-(t-3)/Tu(s)  ds 

This  is  a  concrete  realization  of  the  abstract  problem  (1.2). 
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Equations  of  this  nature  have  been  studied  in  [36]  and  [21],  and  some  numerical 
analysis  is  given  in  [60,  3,  58,  40,  59].  Also,  a  priori  and  a  posteriori  error  estimates  for 
a  finite  element  discretization  of  a  scalar  prototype  ODE  with  memory,  of  the  type  that 
arises  after  spatial  finite  element  semi-discretization  of  this  problem,  are  provided  in  [46]. 


Chapter  5 


Model  problems  in  viscoelasticity: 
internal  variable  formulations 


5.1  Overview  and  notation 

In  this  chapter  we  look  at  how  to  formulate  viscoelasticity  problems  described  by  internal 
variables.  Once  again  our  main  reference  is  Johnson  and  Tessler,  [24].  Before  starting  we 
introduce  some  notation.  For  problems  posed  in  Rn  we  set, 

Hm(Q.)  ~  Hm{n)  x  •  ■  x  Hm(ty  (n  times), 

and,  bearing  in  mind  the  notation  introduced  earlier  in  equations  (4.1)  and  (4.3),  we  set, 

H—[ve  H'iSl)  :  v  =  0  on  Tp). 

Also,  thinking  in  terms  of  weak  formulations  we  also  set, 

(<r,e):=  /  aij£ij(v)  dtt, 

Jn 

and, 

L(t ;  v)  :=  f  v  •  f(t)  dfl  +  <f  v  •  g(t)  dT. 

J  n  J  r  tv 

We  begin  with  quasistatic  problems. 


5.2  Viscostatics:  weak  formulations 


We  recall  first  the  governing  (quasi-equilibrium)  equations  for  quasistatic  viscoelasticity 
introduced  previously  in  (4.3).  Then,  using  the  fundamental  Green’s  formula, 


—  w  jV  dfl  =  /  wv  jdQ,  —  <p  wvrij  dT , 

Jn  J  Jn  J  Jan  3 
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and  also  the  symmetry  of  the  stress  tensor  to  see  that, 

we  then  have  for  v  E  H  that, 

I  fiV{  df2  — —  /  GijjVi  dOi  =  I  (JijVij  dO  (p  (JijTljVi  dT. 

Jo.  Jo,  J  Q,  JrN 

This  implies, 

f  dti  =  f  v  •  /(t)  d£l+  <f  v  •  <j(t)  dl\ 

Jr/vr 

Recalling  the  definitions  made  earlier  we  can  write  this  more  compactly  as, 


(<z,  e(v))  =  L(t;  v)  Vv  e  H. 

This  virtual  work  statement  (weak  formulation)  is  the  standard  “jumping  off  point”  for 
finite  element  approximations  (to  be  considered  later). 

There  now  seems  to  be  (at  least)  two  ways  to  proceed.  The  first  is  to  take  the  route 
described  by  Johnson  and  Tessler  in  [24],  we  describe  this  first. 


Route  1 

We  substitute  for  the  stress  using  the  multidimensional  internal  variables  described  earlier, 
where  we  recall  that, 


Nd 

°ij  “  Cijkl£kl{u(t))  +  XI  ^ijkt 

m— 1 

which  gives  the  representation, 

,  N°r 

(ff,  £(l>))  —  /  (v)  dQ,  +  /  ('ijkl  £kl(t)£ij(v)  dQ- 

Jil  m=l ,/n 


We  now  introduce  internal  displacements  *um  such  that, 


*e%(  *um 


dxi 


+ 


dxk  )' 


and  this  gives  the  virtual  work  statement  in  the  form, 


Nd 

a(u(t),v)  =  L(t;v)  -  £  *am(  *um(t),v) 

m—  1 


Vu  G  H. 
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Here, 

a(u(t),v)  :=  /  Cijki£kl(u(t))eij(v) dll, 

Jft 

*am(  v)  :=  jf  *e£(  *um(*))£i»  dO. 

The  internal  strains  satisfy, 

a  *gg  |  _  ** 

dt  Tm  dt 

and  it  is  important  to  realise  that  these  hold  at  a  point  in  space  and  so  are  effectively 
ordinary  differential  equations. 

Thinking  in  terms  of  the  reversed  internal  variables, 

_ -  *cm 

• —  £  fc  , 

we  may  replace  these  ordinary  differential  equations  by, 

d*eki ,  *eki  ... 
dt  Tm  Tm 

and  solve  these  instead.  This  may  make  more  sense  because  we  know  £  but  not  de/dt. 
Thus  we  arrive  at  a  statement  of  the  problem. 


We  now  describe  the  second  route. 

Route  2 

On  the  other  hand  we  note  that  there  is  no  need  to  introduce  the  internal  displacements 
*um  at  all.  If  we  recall  that, 

nd 

<Jij  =  Cijkl^kl  (u)  d"  'y  ]  &ij  5 
m=  1 
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where, 


* 


'lj 


then  the  virtual  work  statement  takes  the  form, 


nd 

a(u(t),  v)  =  L(t;  v)  -  £  (  *<zm(t),  e(v))  Vv  E  H. 

m=  1 

This  generates  a  different  problem  statement. 

Find  u  :  J  -4  H  such  that, 

ND 

a(u(t),v)  =  L{t\v)  -  £  ( *c rm(t), 6(«))  Vv  €  H. 

TJl—  1 

Here  the  internal  stresses  *<zm  are  given  by, 
d*o*  *o*  _  ^mdeki 

~dT  +  ^T~^kl  dt  ’ 

or:  by  *o*  =  C*$  *e%  where, 

^£M  +  l£M  =  ^i  or  ^£m  +  i£m  =  £M 

dt  rm  dt  dt  Tm  rm 

with  *e$  =  £ki  -  *£$•  _ 


5.3  Viscodynamics:  weak  formulations 


Given  the  quasistatic  problem  it  is  now  straightforward  to  pose  the  viscodynamic  problem. 
Hence, 


nd 

(utt(t),  v )  +  a(u{t),  v)  =  L(t ;  v)  -  ]T) 

m= 1 


*uro(  *«m(f),v) 


WveH. 


The  internal  strain  and  displacement  variables  are  exactly  as  given  earlier.  This  research 
project  is  not  concerned  with  dynamic  problems  and  so  we  leave  the  development  of  this 
area  to  another  time. 
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Part  II 

Numerical  Algorithms  for  stress 

analysis 
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Chapter  6 


Quasistatic  internal  variable 
problems 

6.1  Overview  and  notation 

In  this  chapter  we  consider  some  of  the  possibilities  for  forming  numerical  algorithms 
for  the  quasistatic  problem  modelled  by  internal  variables.  We  consider  semi-  and  fully 
discrete  schemes,  in  that  order,  and  detail  the  computational  implementation.  We  follow 
here  Route  2  as  detailed  earlier  in  Chapter  5. 

First  we  establish  our  notation.  If  W  is  a  function  space  then  for  a  vector-valued 
function  w  :=  (u>j)”=1,  for  which  G  W  for  each  i,  we  shall  write  w  eW  where, 

W  :=  W  x  •  •  •  x  W  ( n  times). 

Moreover,  if  w  :=  (wy)"j=i  is  a  second-order  tensor-valued  function  for  which  u)y  G  W 
for  each  pair  then  we  set, 

\V  =  W  x  •  •  •  x  W  (n  times), 

n  times  n  times 

„ - * - >  , - * - . 

:=  x  •  •  •  x  W )  x  •  •  •  x  (w  x  •  •  •  x  w)  (n  times), 
and  write  w  G  W- 


6.2  Semidiscrete  finite  element  approximation 

We  consider  the  Route  2  problem  in  the  following  form.  Find  u  G  LO0(j7;ff)  and,  for 
each  m,  *am  G  L2(^l))  such  that, 


nd 

a(u(t),v)  =  L(t;v)  -  £  (  Vm(t),e(v))  Vu  G  H, 

m—  1 
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and,  for  each  m, 


d  *am  *em 

dt  +  rm 


Vw  e  L2(Q), 


where  we  define, 


3?  ==  <$?*«(«)■ 

Alternatively,  we  might  decide  to  work  instead  with  the  internal  strains  and  then  this  last 
set  of  constitutive  differential  equations  is  replaced  by, 


(d*em 

dt 


Mw  e  L2  (fi). 


In  this  latter  case  of  course  we  seek  not  only  the  *am  but  the  *em  G  W^(J;  L2(tl))  also. 
(Note  also  that  we  need  u  to  be  time  differentiable,  although  this  requirement  could  be 
removed  by  employing  the  reversed  internal  variables  introduced  earlier  in  Chapter  2) 

To  form  a  notional  semidiscrete  approximation  to  this  problem  we  firstly  partition 
J  :=  [0,T]  into  time  intervals  { Ji }£Ll5  each  given  by, 


Ji  :=  (ti- uU)  with  ki  U  -  U-i 

denoting  the  time  steps.  Now,  for  each  time  interval  Ji  we  partition  Q.  (in  the  usual  way) 
into  a  family  of  disjoint  triangles  suitable  for  piecewise  linear  finite  element  approximation, 
and  then  let  U  denote  the  piecewise  linear  finite  element  approximation  of  u.  We  denote 
the  space  of  piecewise  linear  functions  with  respect  to  this  mesh  (i.e.  during  Ji)  as  Hi, 
and  we  let  flj  represent  the  jth  element  (i.e.  triangle)  in  this  partition. 

In  a  similar  way  we  let  *£m  denote  the  piecewise  constant  (on  each  element)  approx¬ 
imation  to  *am,  and  denote  the  space  of  piecewise  constant  functions  with  respect  to  this 
mesh  (i.e.  during  Ji)  as  Li. 

In  terms  of  the  internal  strains  the  resulting  semidiscrete  finite  element  approximation 
to  this  Route  2  problem  is:  find  U  €  Loo(J-,H ),  satisfying  U\j{  G  L^J^Hi),  and  for 
each  m,  *£m  G  W^iJ;  L2(^l)),  satisfying  *qm\ j.  G  W^JaLi),  such  that  for  each 
i  =  1, 2, . . . ,  N  in  turn, 

nd 

a(U(t),v)  =  L(t;v)  -  £  (  \m(t),e( v))  Vv  6  Hu 

m=  1 


where  *<;Jj  :=  C *e£],  and  the  *e JJJ  are  the  finite  element  approximations  of 
satisfying  (in  each  Ji), 


d*e« 


+ 


-■>w 


dt 


Vw  e  Li . 
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6,3  Semidiscrete  system  equations 


To  represent  a  tensor- valued  function  in  Li  (for  times  t  £  Ji)  we  let  {(frij}  be  a  basis 
consisting  of  piecewise  constant  functions  in  the  following  way, 

_  f  1,  for  x  G  ftj  —  jih  element  in  the  mesh, 
lJ  '  \  0,  elsewhere. 

Then,  faj  €  L{  for  each  j  and  for  an  arbitrary  6  £  Li  we  have  that, 


®kl  —  ^  aj  faj  ? 
i 

where  ay  is  the  (constant)  value  of  the  component  6ki  on  the  element  ftj. 
In  the  expression 


(6,w)  for  w  E  Li, 


we  now  choose  w  such  that  wki  =  <t>ij  is  the  only  non-zero  component,  then, 

(6,w)  =  /  0kiwki  dtl  =  /  dki&jdD, 

Jn  Jfij 

=  meas(fiy)  ay. 

Picking  @  :=  *€m  we  therefore  find  that, 

i_(  •£™,u,)  =  52^, 


and 


meas(fiy)  ay, 


Here  of  course  ay  is  the  spatially  constant  value  of  on  Oy  and  dy 


Also,  since  e(U)  £  we  let  §j(U)  denote  the  spatially  constant  value  of  £(£7)  on 
f lj  and  then  we  may  write, 


3 


Then,  with  the  same  choice  for  w  we  have, 


meas(f&y) 


dpj(U) 

dt 
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Thus,  the  finite  element  approximation  of  the  internal  variable  constitutive  equations 
become, 

dotj  o.j  d/3j(U) 
dt  rm  dt  5 

where, 


otj  is  the  spatially  constant  value  of  *e ^  on  each  fij,  for  every  pair  (fc,Z)  and  for  each  m; 
Pj  is  the  spatially  constant  value  of  Ski{U)  on  flj ,  for  each  pair  (kj). 

To  obtain  the  system  representation  of  the  equilibrium  equations  we  note  first  of  all 

that, 

a{U(t),v)  — ►  AU, 

L(t-,v)  — >  F (t), 

in  a  straightforward  way,  where  A  is  the  elastic  stiffness  matrix  and  F  is  the  load  vector. 
It  now  remains  to  consider  the  viscous  term. 

At  this  point  A.R.  Johnson  would  introduce  the  *Um  internal  displacement  variables 
and  use  the  “Route  1”  formulation  with  the  evolution  equations  written  pointwise  in 
terms  of  the  *Um.  Strictly,  one  should  integrate  these  equations  over  f2  since  derivatives 
of  *Um  are  involved — this  seems  to  introduce  unwanted  complications.  We  consider 
constructing  the  viscous  terms  directly. 

Let  {ipujiy o  be  the  canonical  piecewise  linear  basis  with  respect  to  the  mesh  during 
J7i,  and  let  {'0/}j> o  be  the  resulting  basis  for  H{.  In  the  viscous  term 

(*Sm(t),  6(v)), 


we  choose  v  =  for  some  Z  and  let  S  =  supp^),  then: 

(*sm(f),e(*))  =  E  L 

QiCSJni 

-  Y,  Cijki  *€Slniey(V’()|nimeas(^,), 

SliCS 

=:  F i(t). 

Defining  the  vector  of  viscous  forces  *Fm(t)  :=  (F/(^))^>0  in  this  way  the  system  equations 
take  the  form, 


Nd 

AU (t)  =  F (t)  -  Y  *F m(t), 

m=  1 


daj  aj  _  d(3j(U) 

dt  T-m  dt 


for  all  aj 


on  each  ftj. 


Note  that  *Fm  is  in  principle  easier  to  calculate  than  the  stiffness  contribution  *Am  *Um 
described  by  Johnson  and  Tessler  in  [24],  and  requires  less  memory.  However,  direct  access 
to  the  isoparametric  “builder”  routines  in  the  finite  element  code  is  required. 
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6.4  Equivalence  of  semidiscrete  formulations 

We  now  aim  to  show  that  the  semidiscrete  internal  variable  formulation  described  above  is 
completely  equivalent  (it  has  the  same  solution)  to  the  semidiscrete  Volterra  formulation 
that  we  have  used  previously  (in  for  example,  [45],  [44],  [54]). 

Firstly,  solving  the  ODE  for  aj  in  terms  of  (3j{U)  and  recalling  the  definitions  of  these 
terms  we  have  that, 

**m  =  ek,(U(t))  -  —  f  e-^'T-ekl{U{s))ds 
Tm  JO 

(except  on  element  boundaries),  which  is  expected.  Thus  it  still  makes  sense  to  define, 

:=  ekl(U(t))  -  f  e-^-^euiUWds, 

Tm  JO 

and  we  still  have, 

dt  Tm 

Note  that  by  choosing  w  €  Li  such  that  W{j  =  C^eki{v)  for  any  v  e  Hi  we  have, 

= 


where  we  are  now  working  in  terms  of  the  reversed  internal  strains.  Combining  equations 
then  gives, 


[  *£ijWij  dll  =  f  *e$C$ekl(v)dfl, 

f  *<%ekl(v)dn  =  (*sm,e(v)), 

Jn 


Nd  ( 8  em  \ 

a(U(t),v)  =  L{t;v)  -  ^  Tro(  G  Hi ’ 

m~  1  V  / 

and  where  Wij  :=  C^eki(v). 

Now,  bearing  in  mind  the  combined  equation  given  above  we  note  the  following.  For 
Wki  :=  C^€ij(v)  and  v  G  Hu 

= L  c**1  (* £kiim)  e'(t"s)/rmew(c/(s)) da)  £v{v)  dn' 

=  [  C^iUit^ijHdn-l-  f  f  c:^e-^y^£kl(U(s))etJ(v)dnds. 

JQ  • m  JO  JU 
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Substituting  this  into  the  “combined  equation”  shown  above  then  gives, 

Nd  , 

a(U(t),v)  =  L(t-,v)  -  £  I  C'^euiUmeij^dO 

m= 1  Jtl 
Nd  1  rt  r 

+  £  -  /  /  Cffle-W^eumsVeijWdnda, 

m- 1  Tm  J0  JQ 

which,  using  the  definition  of  a(  • ,  * )  from  earlier,  is  the  same  as, 


/  f  Cijki  +  £  C*$  )  ekl(U(t))£ij(v)  <m  =  L{t\  v) 

Jn  \  m= 1  / 

+  IoIaTs{P1  Cme~{t~S)/Tm^j  eu{U(8))eij(v)dnd8t 

or,  changing  back  to  the  “Dn  notation  for  the  relaxation  function, 

J  Dijki (0)ejtj (U (t) )eij (v)  dS 2  =  L(t;  v)  +  jT  jf  ^  gf  ~  — e«(P~(a))e0-(«)  dfids. 

This  is  precisely  the  Volterra  formulation  we  were  aiming  at.  (Once  again,  this  is  not 
surprising — but  we  have  to  be  certain.)  More  recently  (in  e.g.  [54])  we  write  this  in  the 
form, 

A(U(t),v)  =  L(t; v)  +  [  B(t,s]U(s),v)ds  \/v  € 

J  o 

and  we  will  return  to  this  later.  Now  we  consider  a  fully  discrete  approximation. 


6.5  Fully  discrete  internal  variables 

In  view  of  the  equivalence  just  demonstrated  it  would  be  a  neat  trick  to  find  a  time  dis¬ 
cretization  of  the  evolution  equations  that  is  equivalent  to  our  “Volterra  approximation”, 
and  then  the  existing  error  estimates  given  in  [54,  55]  can  be  applied.  It  is  doubtful 
whether  this  can  be  done.  However,  the  solutions  given  by  the  two  schemes  should  at 
least  be  close,  and  our  next  task  is  to  attempt  to  find  an  expression  connecting  them. 

Let  7 j  be  the  spatially  constant  value  of  on  each  fiy,  for  every  pair  (fc,  l)  in  turn, 
and  for  each  m  in  turn.  Then  by  exactly  the  same  reasoning  as  used  above  to  obtain  the 
ODE’s  for  the  a^s  we  can  also  arrive  at  the  family, 

<hi  +  Tj_  =  0j(U) 

dt  Tm  Tfn 

(Note  that  the  time  derivative  of  / 3j(U )  is  not  required.)  This  is  still  semidiscrete.  To  dis¬ 
cretize  in  time  we  will  use  piecewise  constant  approximations  to  7 j  and  U  denoted  during 
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each  time  interval  Jq  =  (tq-i,tq)  by  %  and  Uq.  A  simple  finite  element  approximation  of 
the  ODE’s  given  above  then  yields,  for  q  =  1,2,..., 


%  -  7?-i  + 


kqlq  _  kqfij(Uq) 

Tm  Tm 


(in  each  f 1j  etc.), 


subject  to  70  =  0)  =  0. 

Setting  rq  :=  kq/rm  we  have  in  general  that, 

(1  +  rq)  %  =  rq(3j(Uq)  +  7,-1, 


and  so  the  first  few  solutions  are  given  by, 


71 

72 


rifijjUi) 

1  +  n  ’ 

rrfjjUi)  nMlA 

l  +  r2  (1  +  r2)(l  +  ri)’ 


~  rjfrjUs)  r2ft(t/2) 

l  +  r3  (l  +  r3)(l  +  r2) 
It  is  not  hard  to  see  that  in  general, 


rxPjjlh) 

(l+r3)(l+r2)(l+ri)' 


~  _  y'  rpPj(Up) 

7?  —  Z_j  q  > 

^Ild  +  r,) 

s=p 

and  that  this  holds  for  every  %  «  7 j  on  Jq  and  where  7 j  —  *e£|  on  Slj.  (Note  that  this 
assumes  the  space  meshes  are  constant  in  time.)  In  tensor  form  we  can  write  this  as, 


V~  1 


rpQjiUp) 

n  (!  +  «-.) 


on  every  element  fij, 


and  where  (*ej%  denotes  the  fully  discrete  approximation  to  ^ej1  during  times  in  Jq. 
Dropping  the  subscript  j  we  therefore  have  that  everywhere  on  the  mesh, 

( *em)q  —  ^2  except  at  element  boundaries. 

P=1II(1  +  ^) 

s=p 


The  question  now  arises  as  to  how  we  determine  the  Up.  For  this  we  use  a  fi¬ 
nite  element  discretization  (in  time)  of  the  semidiscrete  equilibrium  equations.  Using  the 
“combined  equation”  given  earlier  we  then  have,  for  all  suitable  space-time  dependent  test 
functions  v  that, 
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where  we  used  the  semidiscrete  equation  to  write,  for  toy  — 


W  =(e(U(t))~  *em(t),w) 


<9*em 

at  ’ 


and  then  used  the  finite  element  approximations. 

Choosing  the  test  function  v  such  that  it  is  non-zero  only  in  Jq  then  results  in, 

_  I  ftq  No  __ 

o(Ug,  v)  =  -  L(t;  v)dt-  £  (£(Uq)  -  ( *em)q,  w) 

Kq  Jtq-i  m=l 

for  all  v  €  Hq  and  where  toy  =  C*^eki{v).  In  this  case  we  can  write, 

(£(Uq)-(.em)q,w)=  j  C*$ekl{Uq)ei:j{v)<m-  f  C$J(*eS)ffey  («)<«, 

i/ji  j  ri 

and  so  the  fully  discrete  solutions  satisfy, 

[  (  Ciikl  +  7:  Cjjkl  )  eki(Ug)eij(v)  d®1  =  7 ~  f  L(t]v)dt 

\  m= 1  /  ^ 

Nd 

+  £  / 

Using  our  explicit  solution  for  (*em)9  we  then  get, 

f  (  Cjikl  +  yi  Cjhd  )  £kl(Uq)€jj(v)d£l  =  —  f  L(t\v)dt 

Jn  y  m=x  J  K<i  Jk-i 

m=1  P=1II(1  +  ^) 


1  /*^7 

ri 


dt 


Nd  q 

+SSvni=,(i+r.) 


fn  Cij^eki(Up)eij(v)  dQ,. 


Tidying  this  up,  we  find  that  the  fully  discrete  displacements  satisfy, 
f  Dijki(0)£fd(Uq)£ij(v)  —  -r~  (  L{t\v)dt 

Jn  Kq  Jtq- 1 


+ 


5l(n 


Now  we  effect  a  comparison  between  this  fully  discrete  solution  and  the  one  produced 
by  the  “Volterra  formulation”. 
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6.6  Comparison  of  fully  discrete  formulations 


Setting  ip  :=  min{£,£p}  we  can  discretize  the  semidiscrete  Volterra  problem  using  a  piece- 
wise  constant  temporal  approximation  to  the  displacements  to  get, 

A(Uq,v)  =  ^-f  L(t;v)dt  +  -j-  f  YL  [  B(t,s;Up,v)dsdt  Vv  E  Hq. 

kq  Jtq— i  fcq  Jtq—  i  p— i  Jtp— i 

For  the  term  on  the  far  right  we  have, 

(/naA"f  dsdt 

1  Q  rf  /  r  ND  n*'171  \ 

4  'I  '  /  E  dsdt, 

Kq  Jtq- 1  p=1  Jtp-i  \J(l  m=l  Tm  J 

=  -  (  f  f  e  s^Tm  dsdt)  f  Ci]klekl(Up)etj(v) dQ,. 

m= iP=iTmKq  ytq-iJtp-i  /  Jn 

Now,  for  p  =  q  we  have, 

[tq  f  e-(<-s)/Tm  dsdt  =  Tmkq  +  4(e-^/r">  -  l), 

Jtq-lJtq- 1  V  7 


t-g-l  •/  lq 

while  for  p  <  q: 

rta  ft 


rtq  rtp  e~(t-s)/Tm  dsdt  _  T^e -(tq-tp-l)/rm  L kq/rm  _  fg*pAm  _  ^ 

Jtq-!  Jtp- 1  V  y  V  ^ 

Using  these  we  then  have  that  the  fully  discrete  approximation  from  the  Volterra  formu¬ 
lation  satisfies, 

A(Uq,v)  =  ±-  L(t;  v)  dt  +  £  1  + - - -  /  C*$ekl{U  Je^v)  dll 

Kq  Jtq- 1  rn= 1  V  rq  J  Jn 


^  9-1  /e-(tq-tp- 

+  Y  Y  (  ~ 

m~lp=l  \  Q 


—  (ek*/Tm  ~  l)  (efep/Tm  -  l))  Ja  C*$ekl(Up)£ij(v)  dQ. 


Using, 

g(tp-l—tq)/rm  __ 


erp  erp+1  •  •  •  er<*  ’ 
and  expanding  to  first  order  gives, 

e—(tq~tp-l)/Tm 


(ekq/rm  _  _  i)  =  +  +  0(rprq). 


A  similar  approximation  leads  also  to, 

e-fc,/rm  _  j  1  ! 


1  + 


■  +  —  =  T^-+0(r2q). 

rqeTi  1  +  rq  1 
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Thus,  the  fully  discrete  approximation  from  the  Volterra  formulation  also  satisfies, 


J^Dijki(0)eki(Uq)eij(v)  dll  —  ^  L(t',v)dt 

C*ijkCkl(U p)£iJ(v)  dfl 


Nd  q 

+  S|iVn5,(l  +  r .))M' 

Nd  q 

+ 

m=  1  p- 


nd  q  r  _ 

E  EB(rq,rp)  /  C$eu(Up)ev(v)dn, 

m=l  p—\ 


where  E(rq ,  rp)  =  0(rqrp)  is  the  truncation  error  resulting  from  the  approximation  of  the 
exponentials. 


We  can  now  obtain  an  “error  equation”  connecting  these  two  numerical  solutions  as 
follows, 


A{Uq-Uq,v) 


Note  also  that  by  regarding  the  integrals  as  approximations  of  the  rational  functions  we 
could  easily  replace  Up  with  —  Up  in  the  last  term,  and  make  an  adjustment  to  the  first 
term  on  the  right.  That  is,  writing, 


— T “  f  [  e  ^  S^Tm  d S ^  rqrP )? 

Jtq- 1  Jtp- 1 


we  have  also, 


A(Uq-Ug,v)  = 
Nd  q 


EEfrVf  f  e-(t~s)/rm  dsdt\  f  c*j™£kl(Up-Up)etj(v)dn 

m~lp=l  \T™KQ  Jk-1  JtP- 1  /  Jn 

Nd  Q  r  ^ 

-  E  E«(^rp)  / 

m=lp=l 


The  point  about  these  results  is  that  we  can  now  measure  the  distance  between  the  two 
different  numerical  solutions.  This  opens  up  the  possibility  of  using  our  error  estimates 
for  the  Volterra  formulation  to  construct  an  adaptive  solver  for  the  internal  variable  for¬ 
mulations  already  used  by  Dr.  Johnson.  However,  in  the  next  chapter  we  take  a  simpler 
and  more  direct  route  toward  an  adaptive  solver  built  on  the  concept  of  internal  stress 
variables. 


Chapter  7 


The  numerical  algorithm 


7.1  Introduction 


In  this  chapter  we  fix  on  the  numerical  scheme  proposed  in  [54]  for  which  both  a  priori  and 
a  posteriori  error  estimates  have  been  derived  there  and  in  [55,  56].  This  numerical  scheme 
is  based  on  the  Volterra  formulation  of  the  problem  and  discretizes  the  Volterra  integral 
directly.  The  result  is  that  the  entire  solution  history  plays  a  role  in  the  a  posteriori  error 
estimate. 

With  this  in  mind,  below  we  pay  particular  attention  to  marrying  the  theoretical  error 
estimates  derived  for  the  Volterra  formulation  with  a  practical  implementation  drawing 
on  internal  variables.  We  view  this  as  especially  important  since  the  error  bounds  contain 
a  potentially  troublesome  (and  expensive)  term  (denoted  in  the  next  chapter  by  £y).  This 
term  arises  directly  from  the  residual  of  the  finite  element  solution  associated  with  the 
Volterra  integral  on  de-refined  meshes.  It  seems  inconceivable  that  the  effect  of  this 
term  can  be  anything  other  than  to  degrade  the  quality  of  the  error  estimate,  and  it  seems 
possible  that  an  internal  variable  formulation  could  be  used  to  remove  it  completely  with 
comparatively  minor  additional  computational  expense. 

We  say  more  about  the  a  posteriori  error  estimate  in  the  next  chapter  where  we  also 
show  some  adaptive  solutions.  Here  we  outline  the  finite  element  discretization  of  the 
problem  and  quote  from  [54]  the  a  priori  error  estimate,  and  then  give  an  alternative 
numerical  implementation  to  that  in  [54]  by  invoking  internal  stress  variables.  We  term 
the  resulting  scheme  a  hybrid  Volterra-internal- variable  formulation  since  the  numerical 
scheme  and  error  analyses  in  [54,  55,  56]  are  all  framed  in  the  context  of  the  Volterra 
formulation,  while  the  material  below  draws  on  the  internal  variables  discussed  in  the 
previous  chapters  and  is  ultimately  based  on  the  work  of  Johnson  et  al.  in  [24,  25,  23]. 

After  outlining  the  solution  algorithm  we  design  some  test  problems  for  which  we 
know  the  exact  solutions  and  then  compute  some  numerical  solutions,  and  their  errors,  to 
demonstrate  consistency,  spatial  convergence  and  temporal  convergence. 
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7.2  Finite  element  discretization 


To  discretize  the  problem  (4.3)  with  (2.2)  using  the  finite  element  method  we  first  identify 
a  suitable  test  space  H  to  capture  the  spatial  regularity,  and  then  take  the  scalar  product 
of  (4.3)  with  a  test  function  and  integrate  by  parts  over  the  domain  for  an  arbitrary 
time  t  G  J.  Later  we  then  allow  the  test  function  to  be  time  dependent  and  integrate 
through  time  as  well  to  prepare  for  a  space-time  finite  element  discretization. 

Since  the  displacements  are  vector  valued  we  first  define  for  s  =  0, 1, 2, . . .  the  product 
Hilbert  spaces, 

( Hs(n ),  ( ■ ,  •  )s)  :=  (hs(SI),  ( • ,  •  )H>( n))  x  •  •  •  x  (hs{Q),  ( • ,  •  )H*(n))  (n  times), 

where  the  inner  products  are  given  by, 
n 

(w,v)s  :=  53(t«j,«i)jSr.(n)- 
1=1 


for  all  w,  v  G  HS(Q).  These  spaces  have  the  natural  norms  ||  •  ||4  :=  v/(  ■ ,  •  )s  and,  of 
course,  2/2  (^)  =  Also,  and  as  is  usual  for  time  dependent  problems,  for  a  Banach 

space  ( B ,  ||  •  ||s)  we  define  the  Lp(0,t;B)  norms  by, 


IMIl,p(0,t;B) 


(jflK-)lfc*) 


I 

P 


for  t  G  J,  and  with  the  obvious  “ess  sup”  modification  when  p  =  oo. 

Using  the  essential  boundary  condition  we  now  define  the  (spatial)  test  space, 


ff:=  {«eH1(fi):o  =  0onrD},  (7.1) 

and  (see  e.g.  [54]  for  details)  arrive  at  the  “semi- weak”  problem:  find  u  G  L^J^H)  such 
that, 

A(u(t),  v)  =  L(t]  v)  +  f  B(t,s;u(s),v)ds  Vv  G  iif,  a.e.  in  J.  (7.2) 

Jo 

We  call  this  problem  “semi  weak”  because  later  we  will  integrate  over  J  also  to  obtain  a 
“fully  weak”  formulation.  (These  terms  are  just  labels  and  no  mathematical  significance 
should  be  attached  to  the  “semi”  and  “fully”  qualifiers.)  Here  for  the  triangle, 


r  :=  {(*,  s)  G  J  X  J  :  0  <  s  <  t  G  J}, 

the  bilinear  forms  A  :  H  x  H  — >  R  and  B  :  T  *  H  x  H  R  are  defined  by, 

A(w,v)  :=  /  Dijkl(0)eki(w)eij(v)dn, 

Jn 

B(t,s;w,v)  :=  jT  — ^■ekl(w)eij(t ;)dfi, 


(7.3) 


(7.4) 
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for  all  w,  v  £  H,  and  L  :  J  x  H  R  is  a  time  dependent  linear  form  defined  by, 

L(t;v)  :=  /  v  ■  f(t)dtt+  <f  v-g(t)dT.  (7.5) 

Jn  JrN 

For  full  details  on  this  formulation  along  with  some  (rather  standard)  assumptions  we 
refer  to  [54].  In  particular  we  ensure  that  (H,  >!(-,  •))  is  a  Hilbert  space  with  (elastic) 
energy  norm  ||  •  \\h  :=  VM'i  ')  and  dual  H \  and  then  from  [57]  we  have  the  stability 
estimate, 

IMlLp(0,i;tf)  <  S(t)\\L\\Lp(ojt]Ht).  (7.6) 

In  this  estimate  S  :  J  ->  [0,  oo)  is  a  stability  factor  which  for  isotropic  problems  can  be 
expressed  in  terms  of  the  eigenvalues  of  the  stress  relaxation  tensor  (or  matrix)  Q-  This 
factor  is  of  crucial  importance  since  it  appears  in  the  a  posteriori  error  estimates  and 
governs  the  rate  at  which  the  discretization  errors  can  grow  with  time. 

To  carry  out  a  space-time  discretization  we  need  now  to  allow  the  test  function  to 
be  time  dependent  and  then  integrate  the  weak  form  over  the  time  interval  J.  Thus  we 
arrive  at  the  “fully-weak”  formulation  of  this  problem  as:  find  u  £  L^J-,  H)  such  that, 

a(u,v)  =  l(v)  Vn  €  LX(J;H)9  (7.7) 


where: 


a(u,v)  :=  [  A(u(i),v(t))dt  -  [  [  B(t,s;u(s),v(t))  dsdkt, 

Jo  Jo  Jo 

rT 

l(v)  :=  /  L(t;v(t))dt. 

Jo 


(7.8) 

(7.9) 


This  equation  is  the  starting  point  for  the  space-time  finite  element  discretization 
of  the  problem,  and  to  this  we  now  turn.  We  firstly  partition  J  into  N  time  intervals 
where  Ji  :=  (*»-!,<*),  of  lengths  k{  :=  U  -  U-\  >  0  and  such  that, 


J  =  Jx  U  J2  U  ■  •  •  U  JN, 


so  that  to  —  0  and  t n  =  T.  We  use  k  to  denote  the  piecewise  constant  function  such  that 
k\ji  ~  k{. 

For  each  of  these  Ji  we  construct  on  Q  (in  the  usual  way)  a  triangular /tetrahedral 
space-mesh  of  elements  and  denote  the  domain  with  this  mesh  by  ft;.  Element  j  of 
ft;  will  be  denoted  ft;;  and  we  set, 

hij  ~  diam(ft;j), 

and  use  h  to  denote  the  piecewise  constant  mesh  function  given  by  h]^  :=  hij.  We  also 
use  hi  to  denote  the  mesh  function  at  times  t  £  Ji  given  by  hi\n{j  :=  Ay,  and  use  these 
notations  to  see  that  h\j{  :=  hi.  We  need  to  assume  that  arbitrary  partitions  of  ft  x  J  of 
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this  nature  always  exist  because  in  an  adaptive  solution  the  { Ji }  and  {fij  are  of  course 
not  known  in  advance. 

Having  now  broken  the  prismatic  domain  Qx  J  into  the  laminae  (or  “slabs”)  fi  x  J^ 
indexed  by  the  time  levels  i  G  N(l,iV),  we  define  for  each  fi;  the  semidiscrete  (spatial) 
finite  element  spaces, 

Hi  :=  ju  G  H  D  (<7(ft))n  :  v  is  linear  on  Slij  for  each  j  G  N(l,  M*)  j. 

The  space-time  finite  element  spaces  are  now  given  by  Vr  where, 

Vr  :=  {v  G  Lp(J;H) :  v\Jt  G  Vi  G  N(l,iV)}. 

Here  Pr(t^;  Hi)  is  the  vector  space  of  polynomials  of  degree  at  most  r  defined  on  Ji 
with  coefficients  in  H{.  Note  that  our  approximating  functions  in  Vr  are  continuous  in 
space  but  in  general  discontinuous  at  the  knots  These  discontinuities  allow  the 

space-meshes  to  change  with  time. 

We  define  also  the  set  of  internal  edges  (for  triangular  f2 ^  in  R2),  or  faces  (for  tetra¬ 
hedral  O ij  in  R3)  in  each  f as, 

Ci  :=  C  Pi :  3 j  G  N(l,  Mi)  such  that  I  is  an  edge/face  of  |, 

and  the  set  of  edges,  or  faces,  on  the  Neumann  boundary  Tn  as, 

Ti  —  C  Tn  :  3 j  G  N(1,M;)  such  that  t  is  an  edge/face  of  j. 

We  assume  that  there  are  respectively  Nc{  and  such  edges  in  the  time  interval  J{. 
We  also  define, 

Uij  :=  Cj  \  {Li  H  Lj)  for  1  <  j  <  i  <  N,  (7.10) 

as  the  set  of  internal  edges/faces  belonging  to  Cj  but  not  to  Ci .  Note  that  Hu  ~  0,  and 
that  if  we  control  our  adaptivity  and  allow  only  nested  refinements  such  that  \  C  Hi 
for  i  G  N(2,  N ),  then  Ci  0  Cj  =  Cj  and  Hij  =  0. 

Once  we  choose  a  value  for  r,  which  for  us  will  be  r  =  0  or  r  =  1,  we  form  the  finite 
element  approximation  to  (7.7)  as:  find  U  G  Vr  such  that, 

a(U7v)=l(v)  W  G  Vr,  (7.11) 

and  subtracting  this  from  (7.7)  gives  the  fundamentally  important  Galerkin  “orthogonal¬ 
ity”  relationship: 

a(u  —  U,  v)  =  0  Vv  G  VT.  (T.12) 

This  property,  when  coupled  to  the  strong  data  stability  of  an  associated  dual  backward 
problem,  is  the  basic  building  block  in  the  error  estimation  technique  developed  by  Johnson 
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et  al.  in  for  example  [16].  For  our  problem  it  is,  however,  of  limited  use  since  we  do  not 
have  strong  temporal  stability  of  the  solution. 

Note  that  as  the  foregoing  suggests  the  material  in  [54,  55,  56]  is  applicable  to  both 
two-  and  three-dimensional  problems,  although  here  we  are  concerned  only  with  two- 
dimensions.  Also,  we  take  r  =  0  in  all  that  follows  which  corresponds  to  a  piecewise 
constant  temporal  approximation.  The  numerical  algorithm  corresponding  to  piecewise 
linear  time  discretization  is  summarized  in  [54]. 

In  [55]  we  give  a  detailed  list  of  assumptions  on  the  data  and  in  particular  on  the 
approximation  properties  of  the  spaces  {Hi}  and  VT.  These  again  are  rather  standard  and 
would  be  out  of  place  here.  Note  also  that  we  assume  that  D  is  piecewise  constant  in 
space.  We  could  easily  allow  the  case  where  D  is  piecewise  smooth  in  space  at  the  price 
of  an  extra  term  in  the  error  estimate,  but  the  piecewise  constant  case  is  more  likely  to 
arise  in  practical  problems  since  most  materials  have  piecewise  constant  properties. 

The  a  priori  error  estimate  for  this  problem  is  derived  in  [54];  it  takes  the  following 
form. 


Theorem  1  ( A  priori  Galerkin  energy-error  estimate)  Under  certain  natural  as¬ 
sumptions,  and  for  approximation  in  VT,  for  r  =  0, 1,  the  Galerkin  error  e  :=  u  —  U 
satisfies  the  a  priori  error  estimate, 


llu  _  U\\Loo(J;H)  <  C(T)  ^  hD2u 


LootJ-Mm 


+  n* 


fcr+1 


dr+l 


u 


dtr+1 


Lcc{J\H)j 


where  C{T)  is  a  constant.  This  estimate  holds  for  r  =  1  only  if  each  kq  is  small  enough, 
and  depends  also  on  the  ratios  kq/kq-\. 


Such  a  result  is  reassuring  in  that  it  guarantees  convergence  and  also  demonstrates 
the  form  required  for  the  upper  bound  on  the  a  posteriori  estimate  in  order  to  guarantee 
robustness.  However,  from  a  software  implementation  standpoint  the  a  posteriori  error 
estimate  is  much  more  relevant  since  it  can  form  the  basis  of  an  adaptive  code.  We  consider 
adaptivity  in  the  next  chapter  and  below  illustrate  the  convergence  of  the  scheme  with  a 
few  example  exact  solutions.  First  we  detail  the  practical  implementation  of  this  scheme. 


7.3  Internal  variable  formulation 

Before  we  get  to  the  numerical  scheme  in  the  next  section  it  is  useful  first  to  re-write  (7.7) 
in  terms  of  internal  stress  variables.  Recall  first  that  the  constitutive  law  is  given  by, 

cr(t)  =  D(0)e(u(t))  -  f  Ds(t  -  s)e(u(s))  ds 
Jo 

(neglecting  x  dependence),  where, 
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Also,  we  consider  only  isotropic  materials  and  so  write, 

where, 

/I  1  0\  /I  0  0  \ 

I\  :=  [  1  1  0  |  and  I  0  1  0  )  , 

Vo  0  0/  Vo  0  1/2  / 

and  A (t),  fi(t)  are  given  by  Prony  series-type  relaxation  functions.  For  the  stress  we  then 
have, 

tr(f)  =  c(t)  -  S(u;t), 


where, 

<;(t)  :=  D(0)e(u(t))  =  instantaneous  elastic  stress, 

S(tt;  t)  [  Ds(t-  s)e(u(s))  ds  =  inherited  viscous  stress. 

Jo 

Setting  Xi  :=  A and  fa  :=  and  defining  the  internal  stress  variables, 

EAi(u;£)  :=  A ik  [  e~li^~s>)£{u(s))ds 

Jo 

:=  fami  [  e_mi(t_5)e(w(s))  ds, 

Jo 

we  have, 

rt  .  N* 

E(u;t)=  /  Ds(t-s)e(M(s))ds  =  ^SA,(w;t)-l-X)Sw(w;f)- 
Jo  i=  1  i= 1 


Note  that  these  internal  variables  satisfy  evolution  equations  of  the  form, 


and  also  that  the  following  “update  recurrences”  apply, 

SAi(u;f)  =  e~li^~T^'Exi{u-,T)+Xili  J  e~‘i(t~s)e(u(s))ds, 


Using  these  new  definitions  with  (7.8)  and  (7.4)  we  get  an  alternative  representation  of 
the  history  as, 

J  B(t,s;u(s),v)ds  =  j  (/  D5(t  —  s)e(tt(s))  ds'j  •  e(u)df2, 

=  f  H{u\t)  •  £:(v)dO. 

Jn 


The  point  to  note  here  is  that  the  history  integral  has  now  vanished  and  been  replaced 
by  a  local  (in  time)  term.  We  will  use  a  similar  approach  below  to  form  the  numerical 
algorithm. 
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7.4  The  numerical  scheme 

To  obtain  a  practical  scheme  from  (7.11)  in  the  case  r  =  0  we  choose  v  such  that  v  =  0 
outside  of  the  time  interval  Jq,  and  then  v  £  Hq  is  a  piecewise  linear  function  of  x  during 
times  t  £  Jq-  Writing  Uq  :=  U\jq  for  the  discrete  solution  restricted  to  this  time  interval 
equation  (7.11)  then  becomes, 

a(Uq,v)  =  l{v )  Vv\jq  £  Hq,  v\j\jq  =  0, 

for  each  q  =  1,2,....  Here, 

l{v)  =  fq  L(t-,v)dt, 

Jtq- 1 

and, 


a(Uq,v)  =  f  "  A0(t-tq-i;Uq,v)dt-  f  “  f"  B(t,s;U(s),v)dsdt, 

Jtq  —  1  Jtq—  1  J  0 


where, 


A0(t  —  tq- j;  U t?)  D{jki{t 


~~  tq—l)£kl{Uq)£ij{v) 


As  above,  we  assume  throughout  an  isotropic  material  so  that  the  tensor  D  can  be  equiv¬ 
alently  thought  of  (in  two  space  dimensions)  as  the  matrix, 


/  A(i)  +  mW  /*W  0  \ 

D=[  n(t)  A  (*)  +  /!(*)  0  - 

\  0  0  fjt(t)/ 2/ 

Now,  taking  v  to  be  each  basis  function  for  Hq  in  turn,  and  imposing  essential  bound¬ 
ary  data  etc.  in  the  usual  way  we  arrive  at  the  equation  system, 

[Q  A (t-  tq- 1)  dtUq=  [g  F (t)  dt  +  “history” . 

Jtq- 1  Jtq- 1 

We  turn  to  the  “history”  contribution  arising  from  the  double  time  integral  below  since 
we  want  to  give  an  internal  variable  interpretation  along  the  same  lines  as  in  the  previous 
section.  Note  that  this  is  different  to  the  discrete  scheme  presented  in  [54]  where  we 
considered  a  direct  discretization  of  the  Volterra  integral.  In  the  equations  above  A  and 
F  are  essentially  the  same  as  the  standard  stiffness  matrix  and  load  vector  that  arise  in 
standard  finite  element  discretizations  of  the  linear  elasticity  problem.  The  only  difference 
is  that  the  Lame  functions  and  loads  are  time  dependent. 

In  the  practical  scheme  we  integrate  these  equations  to  arrive  at  the  problem:  for 
each  q  =  1, 2, ...  in  turn,  find  Uq  £  Hq  such  that, 


AqUq  =  Fq  +  “history”. 
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Now  A  is  precisely  the  stiffness  matrix  for  linear  elasticity  but  built  with  the  modified 
Lame  functions, 


For  the  generic  Prony  series  relaxation  function, 
N* 

X(t)  =  Xo  +  E  Xie~aiti 
i= 1 


we  have, 


where  kq  :=  tq  —  tg-i.  These  A  and  \i  terms  are  then  easily  evaluated  during  equation 
assembly  by  simple  function  calls.  To  evaluate  the  time  integral  of  the  load  vector  we 
apply  the  following  Gauss  rule  to  the  body  forces  f  and  tractions  g, 


j  m{z)  dz  «  w-  (m(£+)  +  m(£_))  +  w+  ( m(r}+ )  +  m( 77_)) , 
where, 


1  ,  -M  .  1  ,  lh  +  2^30  ,  1  / 15  -  2^30 

ra-’  2 ±V — i40 —  and  ’>±:=2±i i40  • 


This  rule  is  exact  for  m(z)  =  z7,  and  again  we  note  that  the  load  vector  calculation  and 
assembly  routines  in  a  standard  code  need  only  be  modified  in  a  trivial  way. 

To  build  the  “history”  term  into  the  equations  recall  that  we  need  to  include  the 
double  time  integral, 


f  tq  ftq—\ 

/  /  B(t,s-,U(s),v)dsdt. 

J  tq—\  JO 

Recalling  our  internal  stress  variables  from  the  previous  section  note  first  that  we  have, 

/**”*  B(t,s;U(s),v)ds  =  £  f  e-^^UxtiUitg-i)  ■  e(v)dSl, 

Jo  J  n 

N»  t 

+  E  /  •  e(v)dtl. 

»=i  •'n 

From  this  it  follows  that  we  may  write  the  double  time  integral  of  the  strain  history  as, 

rtq  rtq-l 


ftq  Clq~\ 

/  /  B(t,  s;U(s),v)  dsdt 

Jtq—  1  Jo 
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Noting  expressions  of  the  form, 


*  e(v)dfi. 


•  e(v)cifi. 


[tq  dt  =  l~\  1  -  e~likq), 

Jtq-i 

the  double  time  integral  simplifies  to  give, 

r  f1*'1  B(t,s-,U(s),v)dsdt  =  £  [  ■  e(v)dQ 

+  Y,  /  -  e(w)da 

i=l  Jn 

For  a  piecewise  linear  (spatial)  finite  element  approximation  £*.([/; ^_i),  £^(*7;  t9_i) 
and  e{v)  (for  v  6  Hq )  are  constant  on  each  element  and  so  these  terms  are  trivially 
integrated  to  determine  the  local  viscous  load  vectors.  Thus  the  solution  algorithm  takes 
the  following  outline  form. 


Outline  algorithm 

Initialize:  'S\i  =  £/**  =  0  for  each  i. 

Do:  for  time  levels  q  —  1, 2, . . ., 

•  Given 

^Ai(^j^-l)  aild  ^^(Djtg-l) 

from  the  previous  time  step  initialize,  for  each  element  flqj,  the  local  component 
viscous  force  vectors, 

Fl  =  SAi(C7;  V-i)  •  «(»)  and  F%  =  SW(17;  *9_i)  ’  «(«)• 

•  From  these  form, 


(and  similarly  for  the  FqJ. ) . 

•  Form  the  global  component  viscous  forces, 

Fl  =  £  Ft 

fiqj  (Z  Qq 


(and  similarly  for  the  FqJ.). 
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•  Assemble  the  global  total  viscous  forces, 

Ax  A* 

-Fvisc  =  ^2  + 

i=l  i=l 

•  Solve  the  global  system, 

AqU q  =  F  q  +  F  vise- 

•  Update  history  data:  UpdateViscousStressesQ. 

next  q 
Stop 

In  all  of  the  numerical  results  subsequently  presented,  the  global  equations  are  solved 
with  a  diagonally  scaled  conjugate  gradient  iteration.  Unless  explicitly  stated  otherwise 
we  invariably  use  a  residual  tolerance  of  ecg  =  10-7  as  a  stopping  criterion. 

After  the  solution  of  the  global  equations  the  call  to  UpdateViscousStresses() 
performs  updates  of  the  following  form, 

^\i{U-,tq)  =  e~likq'S\i{U-,tq-i)  +  Xik  [tg  e-li^~sh{Uq)ds. 

Jtq-i 

The  viscous  stresses  are  now  ready  for  the  equation  set  to  be  solved  at  the  next  time  level. 
Note  that, 

li  [  q  e~li^tq~s^  ds  =  1  -  e~likq: 

Jtq- 1 

and  so  for  piecewise  constant  temporal  approximation  the  update  equation  simplifies  to, 
Sa  t(U;tq)  =  e-lik<'2Xi(U-,tq-1)  +  \i{l-e-l<k<)e(Uq). 

We  now  detail  some  numerical  tests. 


7.5  Numerical  tests 

In  this  section  we  demonstrate  the  convergence  of  the  numerical  algorithm  by  comparing 
the  known  solution  against  artificial  exact  solutions.  Invariably  we  design  the  loads  and 
tractions  so  that  the  exact  displacements  have  the  form, 

u\(x,y,t)  —  T(t)X(x,y)  and  u2(x,y,t)  =T(t)Y(x,y). 

These  forms  make  it  easier  to  determine  the  loads  /  and  g. 


7.5.  NUMERICAL  TESTS 
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In  the  following  examples  we  use  the  data, 

E  =  2.776  090  556  GPa  and  v  =  0.4,  (7.13) 

which  imply  for  plane  stress  that, 


A(0)  =  1.321 947  884  GPa  and  /z(0)  =  1.982  921 826  GPa.  (7.14) 

The  reason  for  these  choices  will  become  clear  later  when  we  look  at  particular  data  for 
Maranyl  Nylon  6.6.  For  the  moment  we  use  the  arbitrary  relaxation  functions, 

A(t)  :=  A(0)(0.3  +  0.2e-‘  +  0.5e-olt),  (7.15) 

H(t)  :=  /r(0)(0.2  +  0.1e-3i  +  0.2e-0'7*  +  0.3e-2t  +  0.2e-0-2tj.  (7.16) 

Also,  since  the  numerical  results  presented  below  use  various  norms  a  word  on  how  these 
norms  are  calculated  is  in  order. 

For  the  energy  and  2/2  (^)  norms  we  use  a  seven  point  Gauss  rule  (exact  for  quintics — 
see  e#g,  [29]),  while  for  the  L^fl)  norm,  which  we  don’t  take  too  seriously,  we  simply  use 
the  nodal  values.  We  denote  the  resulting  approximate  norms  with  “hats”  as  in, 

II  *  IIjt  ~  II  *  ll§5  II  *  llL2(n)  ~  II  *  ll£2(n)  and  II  ’  IUoo(n)  ~  II  *  ll£00(n)* 

For  the  norms  on  time  dependence  we  approximate  the  local  Loo  norms  by  sampling  at 
the  upper  end  of  the  time  interval  as  in, 


\\w\\Loo(Jq->H)  ~  IMI Loo(Jq;H)  IM^)lltf* 


We  then  approximate  the  global  L0 0  norms  by: 

IMUoo(0,tiH)  «  IMI^o,  t;H)  :=  max  :  0  <  9  ^  *}  forte  Jp. 

The  only  exception  to  these  rules  is  for  the  results  given  for  the  L^J;  Loo  (ft))  norm  of 
the  error  where  we  use  the  following  approximation, 

Mlic^wn))  ~  IMloo  :=  :  9  =  1,2,...  }, 

and  where  tq  ~  (tq  +  tq- 1)/2  denotes  the  midpoint  of  the  time  interval. 

Also,  the  tables  of  results  presented  below  and  in  the  next  chapter  are  verbatim 
output  from  the  finite  element  code  described  later  in  Chapter  9.  The  key  to  the  column 
heading  is  given  in  the  table  below. 
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“verbatim  output” 

= 

“meaning” 

k-i 

= 

time  step  k 

N  elem 

= 

Number  of  elements  M 

1  1 u  I  I _E 

= 

Maximum  energy  norm  IMUoo^/f)  (when  known) 

lluhll-E 

= 

Maximum  energy  norm  \\U\\Loo^j.H) 

1 uh I + | e | 

= 

\\u\\l^(j,h)  +  ||«  -  U\\Loo {J,H)  (when  known) 

|uh|+est | e | 

= 

\\U\\L„(J,H)+£n{T-,U) 

1  |  u-uh  |  |  .£ 

_ 

Maximum  energy  error  ||u  —  (when  known) 

est  I  |e|  I  JE 

= 

Estimated  energy  error  £q  ( T ;  U) 

inf | u-uh I 

= 

Not  used  —  ignore  this  column 

1 1 u-uh | | 

= 

Exact  displacement  error  ||ti  —  ^||z/00(jr;X,2(«))  (when  known) 

1 Is-shl | 

= 

Exact  stress  error  ||<r  -  (when  known) 

7.6  consistency  test 


Using  a  domain  of  arbitrary  shape,  shown  in  Figure  7.1,  we  solve  a  problem  for  which  the 
Galerkin  approximation  is  exact  by  setting, 

T(t)  :=  1,  X(x,y)  :=  —0.03a;  and  Y(x,y)  :=  —  0.04?/. 


For  boundary  conditions  we  impose  a  homogeneous  Dirichlet  condition  on  u  along  the 
left-most  vertical  edge,  and  the  same  on  v  along  the  bottom-most  horizontal  edge.  All 
other  edges  have  tractions  imposed. 

Since  the  approximation  is  continuous  piecewise  linear  in  x  and  y  and  discontinuous 
piecewise  constant  in  t  it  follows  that  the  numerical  solution  should  coincide  with  the 
exact  solution  up  to  quadrature  error. 

We  solve  for  a  constant  time  step  of  k  =  0.5  to  the  final  time  T  —  1  for  mesh  widths 
h  =  0.4, 0.2, 0.1, _ The  first  two  of  these  “regular”  meshes  are  shown  in  Figure  7.1. 

The  results  are  shown  in  the  table. 


ki 

N  .l.m 

1 1 u I |_E 

lluhll.E 

Iuh|+|6l 

iuh|+6stU| 

| lu-uhl l_E 

6st  IMI.E 

inf lu-uhl 

I lu-uhl | 

1 l«-*h| | 

5.0000e-01 

10 

2.6306+03 

2.630e+03 

2.6306+03 

2.6306+03 

1.0377056-06 

2.2983526-07 

1.8446-11 

1.0606-11 

5.4886-02 

5.00006-01 

45 

2.6306+03 

2.6306+03 

2.6306+03 

2.6306+03 

1.0374256-06 

1.5713066-07 

1.8456-11 

1.0606-11 

5.4886-02 

5.00006-01 

177 

2.6306+03 

2.6306+03 

2.6306+03 

2.6306+03 

1.055355e-06 

1.2054936-07 

1.9706-11 

1.0626-11 

5.5326-02 

5. 0000. -01 

762 

2.6306+03 

2.6306+03 

2.6306+03 

2.6306+03 

1.0708046-06 

8.7894046-08 

1.9466-11 

1.0616-11 

5.5556-02 

5.0000e-01 

3109 

2.6306+03 

2.6306+03 

2.6306+03 

2.6306+03 

1.0945676-06 

6.7040616-08 

1.9406-11 

1.0626-11 

5.6016-02 

From  the  table  we  see  that  the  errors  are  small  and  in  fact  at  machine  round-off  once 
the  magnifying  effect  of  A(0)  and  0)  are  taken  into  account. 

Repeating  the  calculations  but  with  k  =  0.25  instead  of  k  —  0.5  we  get  the  following 
tabulated  results. 


ki 

N  6 1 6  m 

llulLE 

lluhll.E 

I uh 1 + 1 6 | 

|uhl+6St(6| 

I lu-uhl |_E 

6Bt  II.II.E 

inf | u-uh | 

1 lu-uhl  | 

ll«-*h|| 

2.60006-01 

10 

2.6306+03 

2.6306+03 

2.6306+03 

2.6306+03 

5.4115076-09 

1.3376376-09 

9.5516-14 

5.5216-14 

3.1736-04 

2.50006-01 

45 

2.6306+03 

2.6306+03 

2.6306+03 

2.6306+03 

1.9354476-07 

4.9638136-08 

1.1596-12 

2.7796-13 

8.3846-03 

2.50006-01 

177 

2.6306+03 

2.6306+03 

2.6306+03 

2.6306+03 

2.0499136-07 

3.9182736-08 

1.4666-12 

3.1856-13 

8.0496-03 

2.60006-01 

762 

2.6306+03 

2.6306+03 

2.6306+03 

2.6306+03 

3.6626636-07 

5.5980726-08 

1.6426-12 

5.0396-13 

1.3696-02 

2.60006-01 

3109 

2.6306+03 

2.6306+03 

2.6306+03 

2.6306+03 

4.5561046-07 

5.1523446-08 

1.8166-12 

5.8806-13 

1.7546-02 

The  estimated  error  quantities  shown  in  the  tables  is  related  to  the  a  posteriori  error 
estimate  and  will  be  explained  fully  in  the  next  chapter.  We  include  it  here  in  order  to  be 
able  to  refer  back  later  on. 


7.7.  SPATIAL  CONVERGENCE 
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Figure  7.1:  Regular  meshes  for  h  =  0.4  and  h  =  0.2 


11  nodes,  10  elements 
10  boundary  sides  during 
times  t  G  (0.500000, 1.000000). 


33  nodes,  45  elements 
19  boundary  sides  during 
times  t  G  (0.500000, 1.000000). 


7.7  Spatial  convergence 

We  now  keep  T(t)  as  above  and  change  X(x,y)  and  Y(x,y)  to, 

X(x,y)  :=  -0.03a:7  +  0.025  sin(27ra:)sin(7ry  +  7r/2), 

Y  (x,  y)  :=  — 0.04a/9  +  0.035  sin(7ra:  +  7r/2)  sm(2ny). 

All  other  data  is  as  before  (with  k  =  0.5)  and  we  again  loop  through  the  same  sequence 
of  regular  meshes  to  examine  the  spatial  error  convergence. 


ki 

N  «l«m 

llull.E 

UuhILE 

|uh|4|*| 

|uh|4«gtl«l 

1 lu-uh 1 I.E 

•ft  i  1*1 I.E 

inf lu-uh I 

I lu-uh I | 

IU-fhll 

6.0000«-01 

10 

8.651*403 

5.856*403 

8.561*403 

9.183*403 

6.244285*403 

7.072830*403 

4.919*-02 

2.393*-02 

2.967*408 

5.0000«-01 

45 

8.561*403 

8.066*403 

8.670*403 

8.689*403 

2.898224*403 

3.233216*403 

1 . 325* -02 

3.972*-03 

1.314*408 

5.0000* -01 

177 

8.561*403 

8.440*403 

8.664*403 

8.594*403 

1.448311*403 

1.616137*403 

7.472«-03 

2.378*-03 

6.703*407 

6. 0000* -01 

762 

8.561*403 

8.631*403 

8.561*403 

8.667*403 

7.154630*402 

7.818353*402 

1.724«-03 

5.423*-04 

3.299*407 

6. 0000* -01 

3109 

8.561*403 

6.554*403 

8.561*403 

8.562*403 

3.422946*402 

3.783819*402 

2 . 925*-04 

5.892«-05 

1.559*407 

5.0000«-01 

12628 

8.661*403 

8.559*403 

8.561*403 

8.561*403 

1.663714*402 

1.856497*402 

5.813*-05 

1.400*-05 

7.579*406 

It  is  evident  that  the  energy  error  is  0(h)  as  expected. 


7.8  Temporal  convergence 

To  examine  the  time  discretization  error  in  isolation  we  now  set, 

T(t)  :=  1  +  i/2  +  sin(27rt),  X(x,y)  :=  -0.03a:  and  Y(x,y)  :=  -OMy. 


We  should  expect  to  see  errors  converge  to  zero  as  k  — ►  0  independently  of  h.  Here  we 
choose  h  —  0.4  corresponding  to  the  mesh  shown  on  the  left  of  Figure  7.1. 


ki 

H  *l*m 

llulLE 

lluhtl.E 

|uh|4 |*| 

|uh|4*gt|*| 

1 lu-uh 1 I.E 

•ft  11*1 I.E 

inf  fu-uhl 

1 |u-uh I  I 

1 Jf-fhl I 

1.0000*400 

10 

3.945*403 

3.178*403 

3.269*403 

3.193*403 

7.677814*402 

3.119966*402 

1.834*-03 

7 . 469*-03 

6.821*407 

6.0000*-01 

10 

3.945*403 

4.641*403 

4.834*403 

4.656*403 

1.986738*403 

7.177276*402 

1.482*-02 

1 . 926*-02 

1.340*408 

2.5000*-01 

10 

6.589*403 

4.764*403 

4.987*403 

4.762*403 

1.798382*403 

6.619429*402 

3.490*-03 

1 . 743*-02 

1.218*408 

1.2500«-01 

10 

6.589*403 

6.404*403 

6.420*403 

5.406*403 

1.060321*403 

3.812408*402 

1.062*-03 

1.018*-02 

7.120*407 

6. 2500* -02 

10 

6.689*403 

5.562*403 

6.663*403 

6.662*403 

5.476366*402 

1.988857*402 

2. 730* -04 

5.309*-03 

3.714*407 

3. 1250*-02 

10 

5.589*403 

6.693*403 

6.693*403 

6.693*403 

2.771164*402 

1.006664*402 

6.807*-05 

2.687.-03 

1.880*407 

1.6626*-02 

10 

6.697*403 

5.595*403 

6.695*403 

6.695*403 

1.390760*402 

5.062766*401 

1.708*-05 

1.348*-03 

9.435*406 

7.8125*-03 

10 

5.697*403 

6.597*403 

6.597*403 

6.697*403 

6.962842*401 

2.629828*401 

4.276*-06 

6. 760* -04 

4.724*406 

3.9062*-03 

10 

6.598*403 

6.697*403 

6.697*403 

5.697*403 

3.483190*401 

1.265696*401 

1.069*-06 

3.377*-04 

2.363*406 

1.963U-03 

10 

6.698*403 

5.698*403 

6.698*403 

6.698*403 

1.741976*401 

6.329462*400 

2.672*-07 

1.689*-04 

1.182*406 
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This  time  we  see  that  the  energy  error  is  O(fc),  again  as  expected. 


7.9  The  general  case:  D  ^  DT 

For  more  general  anisotropic  problems  we  cannot  assume  that  D(t)  —  DT(t)  unless  t  =  0 
or  t  =  oo.  Hence  A  ^  Ar  and  we  would  have  the  additional  complexity  of  solving  a 
nonsymmetric  system  at  each  time  level.  In  this  case  we  could  use  the  evolution  equa¬ 
tions  for  the  internal  variables  and  couple  these  to  an  elasticity  solver  (with  additional 
viscous  loads)  in  an  iterative  solution  algorithm.  This  is  essentially  the  solution  technique 
described  by  Johnson  and  Tessler  in  [24]  and  is  suited  also  to  constitutively  nonlinear 
problems. 


Chapter  8 


Adaptive  error  control 


8.1  Introduction 

In  this  chapter  we  summarize  the  a  posteriori  error  estimates  developed  recently  by  Shaw 
and  Whiteman,  in  [55,  56],  for  the  Volterra  formulation  of  the  quasistatic  problem  defined 
by  (4.3)  with  (2.2).  These  results  follow  on  from  the  exploratory  work  in  [44,  47,  57,  53]. 
We  then  implement  these  error  bounds  in  the  context  of  an  adaptive  space-time  finite 
element  solver  for  the  linear  problem  with  viscoelasticity  described  by  relaxation  functions 
of  Prony  type.  We  note  here  that  while  the  estimates  themselves  do  not  rely  on  this 
form  of  relaxation  function,  it  is  the  only  convenient  choice  for  numerical  computation 
since  it  leads  to  an  economical  history  storage  strategy  (see  for  example  [45]).  Also,  as 
already  demonstrated,  the  Prony  series  leads  to  a  natural  connection  with  internal  variable 
methods. 

The  plan  for  the  chapter  is  as  follows.  We  first  summarize  the  basic  error  bound  as 
given  in  [55],  and  to  do  this  we  refer  back  to  the  weak  form  of  the  problem  as  well  as  its 
space-time  finite  element  discretization  given  in  the  previous  chapter.  Since  in  many  ways 
the  problem  is  close  to  linear  elasticity  we  first  present  the  spatial  error  control  strategy 
in  this  context,  and  also  give  some  numerical  results  on  error  control  via  adaptive  mesh 
refinement. 

In  [55]  we  give  two  forms  for  the  term  in  the  a  posteriori  error  estimate  that  reflects 
the  time  discretization  error.  In  the  first  the  estimate  is  unstable  (as  h  — >  0)  and  so  is 
useless  for  error  control.  The  second  form  on  the  other  hand  is  robust  but  requires  that  the 
residual  be  measured  in  a  discrete  negative  norm.  This  would  require  a  stiffness  matrix 
inversion  and  is  therefore  likely  to  be  prohibitively  expensive. 

These  difficulties  are  due  to  a  lack  of  strong  temporal  stability  in  the  underlying 
differential  equations  which,  in  this  context,  means  roughly  that  there  are  no  time  deriva¬ 
tives  present  on  the  left  of  (4.3)  (as  would  be  the  case  with,  say,  an  ODE  or  parabolic 
equation).  To  overcome  this  fundamental  limitation  (which  has  hampered  other  work  for 
scalar  pure-time  Volterra  equations  in,  for  example,  [28]  and  [2])  we  are  also  working  on 
deriving  a  posteriori  error  estimates  in  a  weak  norm  in  [56].  This  is  based  on  the  prototype 
work  in  [57],  and  we  expect  this  error  estimate  to  allow  for  robust  temporal  error  control 
through  adaptive  time  stepping  as  well  as  the  adaptive  meshing  we  describe  below. 
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8.2  A  posteriori  error  estimate 

In  [55]  we  give  the  following  basic  a  posteriori  Galerkin-energy  error  estimate:  for  each 
discrete  time  t\ ,  £2?  ■  •  •  >  tpi  •  ■ 

llw  ~  U\\Loo(o,tp;H)  —  S(tp)(£n(tp;  U)  +  £j(tp ;  J7)  +  £y{tp\  ^0)>  (8.1) 

where  £^  £j  and  are  residuals  which  are  computable  in  terms  of  the  data  and  the 
finite  element  solution  U ,  and  S(t)  is  the  stability  factor  introduced  before  in  (7.6). 

In  this  section  and  the  next  we  will  be  concerned  only  with  This  term  contains 
the  spatial  discretization  error  (in  the  case  where  'Hij  =  0  for  all  j  <  i)  and  can  be  used 
to  guide  adaptive  space  mesh  refinement.  It  is  essentially  identical  to  the  residual  derived 
for  linear  elasticity  by  Johnson  and  Hansbo  in  [27],  and  this  is  useful  because  in  the  next 
section  we  may  illustrate  its  interpretation  and  use  in  this  less  crowded  context.  The 
extension  to  viscoelasticity  to  come  later  will  then  be  straightforward. 

The  residual  £j  is  the  one  described  earlier  as  being  unstable  (useless)  or— when 
written  in  a  different  form — prohibitively  expensive  to  implement.  As  described  above,  we 
eventually  hope  to  provide  an  alternative  error  estimate  in  which  £q  and  £y  are  essentially 
the  same,  while  £j  is  stabilized  at  the  expense  of  estimating  the  error  in  a  weaker  norm. 

It  is  the  term  £y  that  causes  the  greatest  difficulty  in  this  estimate.  The  spatial 
residuals  in  £q  are  constructed  by  integrating  the  discrete  solution  over  each  element 
to  arrive  at  a  distributional  divergence  of  the  discrete  stress  (compare  (4.3)).  This 
divergence  comprises  two  parts:  the  smooth  function  inside  the  element  (which  is  zero 
in  our  case  of  piecewise  linear  approximation),  and  the  stress  jumps  across  inter-element 
boundaries.  The  difficulty  arises  because  the  stress  is  history  dependent.  This  means  that 
we  have  to  integrate  by  parts  over  not  just  the  elements  in  the  current  mesh,  but  also  over 
all  elements  in  all  previous  meshes.  The  internal  edges  that  appeared  in  previous  meshes 
but  are  no  loner  present  in  the  current  mesh  (e.g.  due  to  derefinement)  are  therefore  “left 
behind”  when  forming  the  standard  residual  /  +  V  *  ah  (which  constitutes  £n),  and  so  we 
consign  the  stress  jumps  across  these  edges  to  the  term  £y.  In  the  particular  case  where 
only  nested  refinements  are  permitted  (so  that  7 Uj  =  0  for  all  j  <  i)  then  no  edges  are 
left  behind  in  this  way  and  we  have  £y  =  0.  This  is  the  case  in  all  of  our  examples  below. 

To  deal  with  mesh  derefinement  would  appear  to  require  fairly  complex  data  struc¬ 
tures  in  the  computer  code  in  order  to  track  all  these  resulting  previous  edges.  Also,  it  is 
not  likely  that  £y  will  act  in  any  way  other  than  to  degrade  the  quality  of  the  estimate 
since  it  contains  historical  contributions  to  the  current  stress.  These  can  then  only  act  to 
reinforce  one  another  in  the  estimate  when  in  fact  the  residual  could  be  much  smaller  due 
to  cancellation.  This  “loss  of  cancellation”  problem  has  been  noted  by  others  in  the  con¬ 
text  of  stress-jump  residual- type  estimators  and  the  memory  in  the  Volterra  integral  acts 
here  only  to  exacerbate  the  problem.  Our  feeling  at  the  moment  is  that  a  representation 
of  the  algorithm  in  terms  of  internal  variables  could  go  some  way  toward  removing  the  £ y 
residual  since  then  all  hereditary  information  is  automatically  represented  on  the  current 
mesh.  The  price  of  this  is  that  the  error  estimates  will  then  be  restricted  to  viscoelasticity 
problems  for  which  Prony  series  relaxation  functions  are  appropriate.  This  does  not  seem 
to  be  an  unreasonable  restriction. 
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We  now  look  in  more  detail  at  the  residual  term  £n-  This  is  defined  at  each  discrete 
time  level  t\,  fy,  ■  •  • ,  tp, . . .  as: 

{tP\  U)  :=  max  {n^J^/IlL^,;/^))  +  n^lKSIUoo^^tfi))}  >  (8-2) 

where  I  Jo,,  and  II  ^  are  constants  appearing  in  certain  interpolation-error  bounds,  and 
hq  =  hq(x)  is  the  piecewise  constant  mesh  function  for  the  mesh  during  times  in  Jq. 

The  first  term  in  the  estimate  is  straightforward  to  interpret  since  it  involves  only  the 
Li{0)  norm  of  the  body  forces  weighted  with  the  mesh  function.  To  define  the  second 
term  we  need  to  establish  various  other  notation. 

We  use  n(tm)  =  (n{'m'l)”=1  to  denote  the  unit  outward  directed  normal  vector  to  the 
the  boundary  <9Qim  of  Qim,  for  m  G  N(1,M;),  and  for  an  edge/face  f  6  £,  we  use  the 
notation  [7*]*  to  denote  the  jump  in  value  across  the  edge/face  i  of  the  components  of 
any  7  =  (7fc)fc=r  That  is,  for  x  G  i  G  A  and  each  i  G  N(l,  N), 

{7k(x)]e  :=  ±  lim  (jk(x  -  enW)  -jk(x  +  en^))  ,  (8.3) 

where:  is  a  unit  normal  vector  to  the  edge/face  t\  and,  to  avoid  elaborate  notation  we 

use  “±”  to  acknowledge  here  that  the  sign  of  this  jump  quantity  is  of  no  interest  at  all. 

We  denote  surface  (edge/face)  integrals  on  the  element  boundaries  by 
(w,v)e:=  w-vdl  with  II  •  III  ■•=  \Z(-,  ■)*> 
and  we  define  the  discrete  traction, 

mum  :=m-,umu, 

where  for  a  unit  vector  n, 

9k(t]U(t))  ^Dkuj(0)eij(U(t))  —  £i j(U{s))ds^ni  (8.4) 

(compare  the  natural  boundary  condition  in  (4.3)).  In  general  these  discrete  tractions  will 
not  be  uniquely  defined  on  any  edge/face  I  G  A,  but  the  jumps  {gje  will. 

Now,  for  each  time  level  we  define  r\jq  =  (rk)k=i  by, 


\\l9k(^U{t))]t\,  ioxItCq, 
rk(t',  U(t))  :=  <  | gk{t)  -  gk(f,  t/(t))|,  on  TN, 

.  0,  on  r0, 

and  then  with  all  of  these  definitions  we  define  S  £  Li(0)  by, 


SI n„-  := 


l|r(t;t/(t))k2(an„) 

yJhqjmeas(Qqj) 
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for  all  j  6  N(l,  Mq)-  Note  that, 


E 

\nqj  cfi. 


hqjT(t\U  (t)) 


1 

2 


(8.5) 


We  now  use  £q  to  derive  element  error  indicators  and  a  mesh  refinement  criterion  in 
the  context  of  linear  elasticity. 


8.3  Adaptive  meshing  for  linear  elasticity 


In  the  linear  elasticity  context  there  is  of  course  no  time  dependence  and  so  we  get  the 
simpler  result, 

||w  -  U\\H  <  £q{U)  :=  nn||/^f  ||l2(q)  +  Ilt\\hS\\L2(n)- 

Note  that  the  Loo  norms  and  the  subscripts  p  and  q  marking  the  time  levels  are  not 
needed  here.  Also  we  have  temporarily  set  S(t)  =  1.  It  will  be  straightforward  to  extend 
the  following  results  to  the  viscoelasticity  problem  later. 

The  goal  is  to  design  a  mesh  for  which  \\u  —  U\\h  <  TOL  where  TOL  >  0  is  a 
user-defined  tolerance  level.  Further,  this  mesh  should  be  optimal  in  the  sense  that  the 
error  control  is  achieved  with  as  few  degrees  of  freedom  as  possible.  Although  impossibly 
difficult  to  obtain  in  an  exact  sense  (see  for  example  the  discussion  in  [17]),  such  a  mesh 
can  be  approximated  with  sensible  and  restrained  use  of  an  a  posteriori  error  estimate. 
We  aim  for  a  mesh  size  modification  strategy  for  each  element  Qj  in  the  mesh. 

The  error  control  is  clearly  achieved  if  we  ensure  that  £q(U)  <  TOL,  and  to  do  this 
we  use  the  technique  of  equidistribution  whereby  each  element  is  allowed  to  contribute 
equally  to  the  global  error  (regardless  of  element  size).  In  practice  this  means  that  we 
assign  a  local  tolerance  to  I  >  0  to  each  element  and  attempt  to  control  the  local  element 
errors  to  within  tol. 

For  each  element  Qj  in  the  mesh  we  seek  a  local  meshwidth  hj  and  solution  U  for 
which, 

v]  :=  n^ll/lli^.)  +  n^usiil^.)  =  tol. 

Rearranging  this  we  then  arrive  at  an  adaptive  mesh  size  selector  via, 


Lnew  _ _ 

n3  •_ 


tol 


+  n?l|9||2i2(n.) 


(8.6) 


for  each  element  f lj.  Then,  globally,  where  M  is  the  number  of  elements  in  the  mesh  we 
get, 

M  M 

M  x  tol  =  tol  =  E7?!’ 

j=]  j- 1 
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=  nl\\hf\\l2{Q)  +  nt\\h5\\l2{n)  >  ±(«h(to)2, 

where  we  used  the  inequality  a2  +  b2  >  (|a|  +  |6|)2/2.  Hence, 


£n(U)  <  TOL 


is  guaranteed  if 


to  I  := 


TOL2 
2  M  ’ 


This  is  then  our  formula  for  the  local  error  tolerances. 

The  basic  form  of  the  adaptive  algorithm  we  implement  is  then  as  follows. 
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1 .  Set  i  —  0  and  generate  an  initial  mesh  JVC1. 

2.  Solve  for  the  displacements  Ul  on  M\ 

3.  Compute  rf-  on  each  element  tlj  then: 

if  <TOL2/2  STOP 

else  determine  a  new  mesh  Ml+1  from  (8.6) 

4 .  Set  i  =  i  +  1  and  repeat  from  Step  2 . 


At  this  point  we  should  sound  a  note  of  caution.  This  algorithm  should  not  be 
interpreted  glibly  in  that  the  else  clause  in  Step  3  should  not  be  taken  literally.  If  every 
element  is  refined  according  to  (8.6)  then,  unless  the  initial  mesh  M°  is  already  nearly 
optimal,  one  is  likely  to  get  severe  over-refinement  in  the  adapted  meshes  due  to  pollution 
error.  Instead,  rather  than  attempt  to  derive  the  required  mesh  in  a  “one  shot”  manner 
by  subdividing  each  elements  as  indicated  by  (8.6),  one  should  only  halve  the  size  of  the 
elements  marked  for  refinement  and  then  loop  back  and  recompute. 

Before  moving  on  to  viscoelasticity  we  now  give  a  couple  of  examples  to  illustrate  the 
adaptive  solution  of  linear  elasticity  problems.  We  again  assume  an  isotropic  material  in 
which  case  Hooke’s  law  becomes, 


/A  +  /x  A  0  \  /  en  \ 

I  A  A  +  /i  o  I  |  €22  I  > 

Vo  0  /i/2/  \2enJ 

where  A  and  p  are  the  Lame  coefficients  related  to  the  Young’s  modulus  E  and  Poisson’s 
ratio  v  through, 


A  := 


vE 

(1  +  u)(l  -  2v) 


and 


/i  := 


E 

l  +  u 


Note  that  also  \x  —  2 G  where  G  is  the  shear  modulus  of  the  material.  In  two-dimensional 
problems  in  which  the  stress  in  the  third  direction  is  negligible  engineers  frequently  employ 
the  plane  stress  approximation  where, 
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The  first  definition  of  A  in  a  two-dimensional  formulation  implies  plane  strain  (in  which 
case  the  strain  in  the  third  direction  is  assumed  negligible).  Typically  one  may  use  plane 
stress  for  very  thin  components,  and  plane  strain  for  very  thick  ones. 

In  the  following  examples  we  use  the  same  data  as  before:  (7.13),  (7.14),  (7.15)  and 
(7.16). 

8.3.1  The  interpolation-error  constants  Iln  and  11^ 

To  complete  the  description  of  the  a  posteriori  error  estimate  for  linear  elasticity  we  need 
to  specify  the  interpolation-error  constants  Hq  and  II^.  For  this  we  adapt  the  approximate 
values  calculated  by  Ludwig  in  [33,  Tables  6.2  and  6.4].  Here  the  interpolation  constants 
are  given  for  a  variety  of  Poisson  ratios,  for  plane  strain  with  E  =  1,  and  assuming  a  mesh 
of  right-angled  triangles.  The  values  are  reproduced  here  in  Table  8.1. 


Table  8.1:  Interpolation-error  constants 


V 

Iln 

n* 

Iln /Ilf 

0.1 

1.143177 

2.683068 

0.426071 

0.2 

1.172415 

2.715093 

0.431814 

0.3 

1.197708 

2.734093 

0.438064 

0.4 

1.219404 

2.741196 

0.444844 

0.5 -e 

1.237806 

(2.947157) 

0.42* 

The  value  for  for  v  — >  0.5  is  not  given  by  Ludwig.  To  calculate  it  here  we  assume 
a  scaling  of  Hq/He  =  0.42  and  then  work  from  the  tabulated  value  of  Hq.  These  values 
can  be  incorporated  into  a  computer  code  as  a  look-up  table  and  then  values  for  any  value 
of  v  obtained  by  linear  interpolation.  The  constants  scale  with  E  in  the  following  way: 


Hn  _  Lb  __ 

nn  ~~  fit  ~  V  e' 


(Note  also  that  Ludwig  uses  2/x  in  Hooke’s  law  where  we  use  //,  but  since  this  does  not 
affect  E  and  v  this  difference  is  immaterial). 


We  need  now  to  obtain 
note  first  that: 


A  = 
V  = 


vE 


(1  —  2z/)(l  -f  i^) 
E 


1  +  v 


> 


the  corresponding  interpolation  constants  for  plane  stress, 

'  ^  _  {a + 3a)/^ 

=>  <  ^  for  plane  strain, 

^  fi  +  2A 


A  = 


V  = 


uE 


(1  -  ^)2 
E 

1+is 


E  = 


i/  = 


(ju  -b  2A)^i 

- b  A 
A 


for  plane  stress. 


/i  +  A 

So,  given  E  and  v  for  plane  stress  we  can  first  calculate  A  and  /z,  and  then  work  backwards 
with  these  values  to  find  the  corresponding  E  and  v  for  plane  strain — E  and  V  say.  From 
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these  we  can  then  determine  the  interpolation  constants  from  the  table  and  the  scaling 
given  above. 

In  our  numerical  results  we  take  v  —  0.4  and  E  =  2.776  090  556  GPa  which  give  the 
plane  stress  values, 


A  =  0.476 190  476 E  and  y.  =  0.714  285  714 E. 

These  give  the  corresponding  E  and  V  plane  strain  values, 
E  =  0.918  367  3461?  and  V  =  0.2857 .... 


Assuming 


E=l 


the  table  gives  for  this  Poisson  ratio, 


ftn  «  1.2  and  fb  ~  2.8, 


and  then  using  the  scaling  for  E  we  finally  get  the  values, 


_  1.252198 

=  ~^Je 


and 


nt  = 


2.9218 

y/E  ' 


However,  there  is  some  doubt  as  to  whether  this  exercise  is  worthwhile.  These  constants 
represent  the  worse  possible  case  in  interpolation  error  and  often  end  up  making  the  a 
posteriori  error  estimate  significantly  over-estimate  the  finite  element  error.  To  address 
this  difficulty  we  have  calibrated  these  constants  against  exact  solutions  in  order  to 
render  the  a  posteriori  estimates  more  realistic.  The  end  result  is  that  we  divide  the 
values  given  above  by  a  factor  of  ten  and  twenty  respectively.  These  values  are  suggested 
by  the  calculations  for  the  exact  solution  given  in  an  earlier  table,  but  we  will  see  below 
that  it  is  desirable  to  determine  a  more  systematic  calibration  technique. 


8.3.2  Example:  exact  solution 

For  the  moment  we  switch  off  all  time  dependencies  and  viscoelasticity  effects  and  consider 
the  adaptive  solution  of  a  linear  elasticity  problem.  We  use  the  elastic  coefficients  A  :=  A(0) 
and  \x  :=  fi( 0)  from  the  previous  chapter,  and  impose  loads  and  tractions  such  that  the 
exact  solutions  are, 

T(t)  :=  10-2,  X(x,y)  :=  —0.03a;20  and  Y(x,y)  :=  -0.04y20. 

We  need  these  because  adaptivity  doesn’t  really  show  up  as  useful  for  smooth  solutions — so 
we  simulate  a  singularity. 

The  notation  here  is  exactly  as  in  the  previous  chapter  and  we  impose  boundary 
conditions  in  the  same  way.  For  these  numerical  tests  we  simply  solve  to  time  T  =  1  using 
a  single  time  step  k  —  1.  The  first  set  of  results  show  reference  solutions  calculated  for 
uniform  meshes  with  h  =  0.4, 0.2, _ Again,  these  are  exactly  as  in  the  previous  chapter. 
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ki 

N  *l«m 

llull.E 

lluhl I.E 

|uhl+l« 1 

|ub|+**t |*| 

| |u-uh| | _E 

•«t  IUILe 

inf | u-uh | 

I lu-uhl | 

1 l»-*h| | 

1. 0000*400 

10 

8 . 023*+01 

5 . 224«+01 

8. 198*+01 

1.380* +02 

6.317922*+01 

1.277228*402 

6.289*-04 

2. 822* -04 

3.896«+06 

1.0000* +00 

45 

8.370*+01 

6.688«+01 

8.377*+01 

9.146«+01 

5. 173867«+01 

6 .344186* +01 

1.402*-04 

4.603*-05 

3.109*+06 

1.0000®+00 

177 

8.397*+01 

7.671*+01 

8. 397* +01 

8.380*+01 

3.416709*+01 

3. 374743* +01 

1.086«-04 

2. 062* -05 

2.098*406 

1.0000*+00 

762 

8 . 398«+01 

8. 167*+01 

8.398*+01 

8.342*+01 

1.954781*401 

1.700767*401 

2.818«-05 

6. 674* -06 

1.186*406 

1.0000* +00 

3109 

8.398*+01 

8.347*+01 

8.398e+01 

8.386*+01 

9.258688«+00 

8.078055*400 

1.351*-05 

2.936*-06 

6.496*405 

1.0000* +00 

12628 

8 . 398*+01 

8.387«+01 

8.398*+01 

8.396*+01 

4. 291128«+00 

3.861923*400 

2.782«-06 

7. 406* -07 

2.656*405 

Looking  at  the  estimated  errors  we  now  tabulate  adapted  solutions  for  TOL  =  15, 
12.5,  10,  7.5,  5,  with  h  —  0.1  for  the  initial  mesh. 


ki 

N  *l«m 

llull.E 

lluhl I.E 

I uh |+|*| 

|uh|+«st|*| 

1 lu-uhl I.E 

•st  IUILE 

inf | u-uh | 

I lu-uhl | 

1  U-#h|  | 

1.0000«+00 

488 

8.398*401 

8.300*401 

8.398*+01 

8.376*401 

1.275528*401 

1.124042*401 

1.806«-05 

3.435«-Q6 

7.766*405 

1.0000*400 

639 

8.398«+01 

8.325«+01 

8.398«+01 

8.382*401 

1.103850*401 

9.734371*+00 

1 . 838*-05 

2.866*-06 

6.731*405 

1.0000«+00 

711 

8.398*401 

8.333*401 

8.398«+01 

8.386*401 

1.042688*401 

9.379712*400 

1.772*-05 

2. 634* -06 

6.384*405 

1.0000*400 

1085 

8.398*401 

8.360*+01 

8.398*+01 

8.389*401 

7.918571*400 

6.890573*400 

9.779*-06 

1.720*-06 

4.819*405 

1.0000*400 

1858 

8.398*401 

8.377*401 

8,398*+01 

8.393*401 

5.874564*400 

5.137793*400 

6.566*-06 

1.122«-06 

3.582*405 

Repeating  these  calculations  but  starting  with  an  initial  mesh  of  h  =  0.2  gives  the 
following  results. 


ki 

N  *l«m 

llull.E 

lluhll.E 

(uh 1+1*1 

|uh| 4«ct | • | 

| |u-uh| | _E 

•st  IUILE 

inf | u-uh | 

f lu-uhl | 

IU-sh|  | 

1.0000*400 

512 

8.398*401 

8.302*401 

8.398«+01 

8.376*401 

1.263565*401 

1.099203*401 

2.475*-05 

3.259*-06 

7.600*405 

1.0000*400 

561 

8.398*401 

8.311*401 

8.398«+01 

8.375«+01 

1.206217*401 

1.036425*401 

2.377*-05 

2.912*-06 

7.248*+05 

1.0000«+00 

749 

8.398*+01 

8.333*401 

8.398«+01 

8.384*401 

1.045034*401 

9.294968*400 

1 . 896*-05 

1.916*-06 

6.337*405 

1.0000*400 

1289 

8.398*401 

8.361*+01 

8.398*401 

8.388*401 

7.842236*400 

6.718701*400 

1 .051* -05 

1.479* -06 

4.728*405 

1.0000*400 

2403 

8.398*401 

8.379e+01 

8.398*401 

8.393*401 

5.561009*400 

4.808581*400 

6. 485* -06 

7.713*-07 

3.352«+05 

Although  the  effect  on  the  number  of  elements  is  not  too  serious,  this  seems  to  verify 
the  popular  wisdom  that  the  effectiveness  of  the  adaptive  meshing  procedure  is  influenced 
by  the  initial  mesh.  We  will  show  some  pictures  of  adapted  meshes  in  the  next  section 
where  we  move  on  and  perform  the  analogous  numerical  experiments  for  linear  viscoelas¬ 
ticity. 

8.4  Adaptive  meshing  for  linear  viscoelasticity 

With  the  strategy  for  adaptive  space  meshing  established  in  the  context  of  the  linear  elas¬ 
ticity  problem,  it  is  now  straightforward  to  extend  it  to  the  time  dependent  viscoelasticity 
problem.  All  that  we  do  is  apply  the  mesh  refinement  criteria  at  each  time  level  in  turn,  as 
if  we  were  dealing  with  a  sequence  of  linear  elasticity  problems.  The  only  major  difference 
is  the  inclusion  of  the  time  dependence  in  the  a  posteriori  error  estimate.  Recalling  this 
estimate  from  (8.1)  we  now  assume  that  we  are  given  a  tolerance  TOLp  with  which  to 
control  the  mesh  such  that, 

S(tp)£n(tp ;  U)  <  TOLq,  for  each  time  level:  p=  1,2, - 

We  describe  the  stability  factor  in  more  detail  below  and  for  now  note  only  that  it  is  a 
non-decreasing  function.  So,  replacing  S(tp)  with  S(T)  and  recalling  the  definition  in  (8.2) 
we  see  that  this  error  control  will  be  guaranteed  if  we  ensure  that, 


^qWhqf\\boo(Jq\L2(Ci))  +  n/||ft9S|Uoo(Ji;^2(n))  -  ~S{T) 
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at  each  time  level  tq.  The  strategy  is  now  exactly  as  above  for  the  linear  elasticity  problem. 
Assuming  a  local  tolerance  tol  we  derive  the  adaptive  mesh  size  selector, 


i  new _ 

aqj  . 


tol 


ll/ll2 


+  n?||3||2 


(8.7) 


and  then  S(tp)£n(tp,  U)  <  TOLn  is  guaranteed  if  we  choose 

tol£ 

t0l-2. MSP(TY 

Below,  in  the  implementation,  we  replace  the  indicated  norms  with  practical  approxima¬ 
tions  as  discussed  earlier.  Note  that  in  the  above  we  have  set  Ifo,  =  nn  for  each  q:  i.e. 
we  take  the  interpolation-error  constants  as  time  independent  (even  though  the  mesh  is 
not),  and  then  use  the  values  for  Iln  and  II*  given  previously. 

We  now  outline  the  form  of  the  stability  factor  and  then  follow  with  the  extension  of 
the  numerical  results  for  linear  elasticity  to  this  time  dependent  problem.  For  the  moment 
we  assume  that  there  are  no  time  discretization  errors  since  we  want  to  consider  the  more 
difficult  problem  of  temporal  error  control  in  a  separate  section  below. 


8.4.1  The  stability  factor  S(T) 

For  an  isotropic  viscoelastic  material  in  two  dimensions  the  stability  factor  S(t)  has  been 
derived  rather  precisely  by  Shaw  and  Whiteman  in  [57].  We  give  here  only  the  main  result 
and  refer  to  the  reference  for  details. 

For  a  viscoelastic  solid  we  have, 


S{t):=(l-J  <f>(s)ds^J  where  <f>(t)  :=  max{wi(t),o^(f)}, 


and, 


u>i(t)  :=  - 


A  '{t)  +  n'(t) 


W2 (<)  :=  - 


At) 


A(0)  +  n{0)  ’  p(0 )' 

For  example,  using  our  test  data  from  (7.15)  and  (7.16)  we  have, 

_  r  W2(t),  for  0  <  t  <  t*  :=  2.28476. . ., 

|  a;i(f),  for  t  >  t*. 


these  give, 


S(t)  :=  { 


m(Q) 

nit)’ 


for  0  <  t  <  t\ 


U(*(0)  A(0)  +  /i(0)  J 
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Note  also  (thinking  ahead  to  the  “Maranyl”  Nylon  6.6  data  given  later)  that  for  a  syn¬ 
chronous  (solid)  viscoelastic  material  there  exists  a  generic  relaxation  function  ip(t),  nor¬ 
malized  to  (p( 0)  =  1,  and  such  that, 


m 

A(0) 


p{i) 


In  this  case  we  have  the  simpler  result  S(t)  =  1 


8.4.2  Example:  exact  solution 


We  continue  with  the  “singular”  linear  elasticity  example  used  above,  but  now  use  the 
assynchronous  relaxation  data  for  A (<)  and  n(t)  as  described  by  equations  (7.13) — (7.16). 
In  particular  we  still  take  T(t)  =  10-2  and  solve  up  to  time  T  =  1,  this  means  that  the 
solution  is  time  independent  and  the  only  time  discretization  error  is  due  to  numerical 
integration  of  the  load  terms  f  and  g.  To  keep  this  quadrature  error  “under  control”  we 
first  do  some  tests  to  determine  an  appropriate  time  step. 

To  follow  the  pattern  of  the  linear  elasticity  calculations  we  first  show  results  for 
uniform  meshes  h  =  0.4, 0.2, ...  in  the  tables  below  for  k  =  1.0,  0.5  and  0.25.  The  first 
table  is  for  k  =  1.0. 


ki 

K  a  lam 

llull_E 

lluhll.E 

|uh|+ |a| 

|uh|+aat|a| 

| |u-uh| 1 _E 

aat  j | a | | _E 

inf | u-uh | 

Mu-uh|| 

1 1 a-ahl | 

1.0000e+00 

10 

8.023a+01 

5.237a+01 

8.209e+01 

1 . 686a+02 

6.321785e+01 

1.602739a +02 

5.074a-04 

2.774a-04 

2.584a+06 

1.0000a +00 

45 

8. 370a +01 

6.616a+01 

8.400a+01 

1.032a +02 

5.176530a+01 

7.917671a+01 

1.777a -04 

5.827a-05 

2.041e +06 

1 . 0000a+00 

177 

8.397e+01 

7.677a+01 

8.404a+01 

8.760a+01 

3.417289a+01 

4.217796a+01 

1 . 328a-04 

3.048a-05 

1.391a+06 

1.0000a+00 

762 

8.398a+01 

8. 171a+01 

8.401a+01 

8.442a+01 

1.955150a+01 

2.121986a+01 

3.469a-05 

9.753a-06 

7.834a+05 

1 .0000a +00  3109  8.398e+01  8.348a+01  8.399a+01  8.408a+01  9.260549a+00 

The  analogous  results  for  k  —  0.5  now  follow. 

1.005926a+01 

1.368a-05 

2. 587a -06 

3. 615a+05 

ki 

N  alam 

llull.E 

lluhll.E 

luhl+fa 1 

|uhl+ast|a| 

1 1 u-uh 1 1 .E 

a at  llall.E 

inf  1 u-uh 1 

| lu-uhl | 

1 la-ahl 1 

S.0000e-01 

10 

8.023e+0i 

5. 242a +01 

8.214a+01 

1 . 686a+02 

6. 323676a +01 

1.602618e+02 

5. 159a-04 

2.791a-04 

3 . 008a+06 

B.OOOOa-Oi 

45 

8.370a+01 

6.625a+01 

8.408a+01 

1 . 032a+02 

5. 176755a +01 

7.914604a+01 

1.893a-04 

6.411a-05 

2.384e+06 

5.0000a -01 

177 

8.397a+01 

7.680a+01 

8.406a+01 

8.761a+01 

3. 417703a +01 

4.216685a+01 

1.403a-04 

3.398a-05 

1.619e+06 

6.0000a -01 

762 

8. 398a+01 

8. 172a+01 

8. 403a +01 

8. 443a +01 

1.955415a+01 

2.121487a+01 

3.670a-05 

1.085a-05 

9.129a+05 

5 .0000a -01  3109 

And  now 

8 . 398a+01  8.348a+01 

the  results  for 

8.400a+01 

k  =  0.25. 

8. 409a +01 

9.261886a+00 

1.005697a+01 

i.372a-05 

2.724a-06 

4.219a+05 

ki 

N  a  lain 

llull.E 

lluhll.E 

juh|+|e I 

luhl+aat |a 1 

I lu-uhl 1 _E 

aat  ||a||_E 

inf | u-uh | 

1 lu-uhl | 

1 l«-ah| 1 

2.5000a-01 

10 

8.023e+01 

6.244a+01 

8.216e+01 

1.686a+02 

6. 324584a+01 

1.602583a+02 

6. 217a -04 

2.804a-04 

3.361a+06 

2.5000a-01 

45 

8 . 370a+01 

6.629e+01 

8.411a+01 

1.032a+02 

5.177362a+01 

7.913404a+01 

1.942a-04 

6.671a-05 

2.671a+06 

2.5000a-01 

177 

8.397a+01 

7.681a+01 

8.407a+01 

8.762e+01 

3.417912a+01 

4.216233a+01 

1.435a-04 

3.648a-05 

1.809a+06 

2. 5000a-01 

762 

8.398a+01 

8.172a+01 

8.403a+01 

8.443a+01 

1.955547a+01 

2.1212836+01 

3. 754a -05 

1 . 132a-05 

1.021a+06 

2.6000e-01 

3109 

8.398a+01 

8.349a+01 

8.400e+01 

8. 409a +01 

9.262551a+00 

1.005604a+01 

1.375a-05 

2.818a-06 

4.725a+05 

One  can  see  that  the  quadrature  error  has  only  a  marginal  affect  on  the  “energy” 
quantities,  which  are  our  primary  concern,  and  hence  we  use  only  a  single  time  step  for 
the  following  adaptive  solutions. 

For  the  adaptive  solutions  we  use  the  same  set  of  tolerances  as  before:  TOL  =  15, 
12.5,  10,  7.5,  5.  The  initial  mesh  is  again  given  by  h  =  0.1,  and  the  results  are  shown  in 
the  table  below  (with  k  =  1). 
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ki 

N  #l#m 

llulLE 

lluhll.E 

|uh|+|#l 

|uh|+##t|#| 

I lu-uhl | _E 

««t  M.ILE 

inf | u-uh | 

1 fu-uhl | 

1 1 ■-•hi  I 

1.0000#+00 

695 

8. 398# +01 

8.331«+01 

8.398#+01 

8.415#+01 

1.061322#+01 

1. 186403#+01 

2. 015# -05 

3. 930# -06 

4. 303# +05 

1 . 0000e+00 

711 

8 . 398#+0i 

8.333#+01 

8.398#+01 

8.415#+01 

1.042769#+01 

1 . 171360#+01 

1.986#-0S 

3. 928# -06 

4.234«+06 

1.0000#+00 

1056 

8.398«+01 

8.359#+01 

8.398«+01 

8.405#+01 

8.086324#+00 

8. 724684# +00 

1.217#-05 

2.726#-06 

3.261#+05 

1.0000«+00 

1584 

8.398#+01 

8. 373# +01 

8. 398# +01 

8.402#+01 

6. 478367# +00 

7. 022669# +00 

8.442#-06 

2.027#-06 

2.612e+05 

1.0000# +00 

3064 

8.398#+01 

6.385«+01 

8.398#+01 

8.400«+01 

4. 699467# +00 

5.016549«+00 

4.422#-06 

1.241#-06 

1.884#+05 

The  adapted  meshes  for  a  selection  of  tolerances  and  some  plots  of  the  stresses  are 
shown  in  the  Figure  8.1. 

This  example  is  of  course  non-physical  and  is  useful  only  because  we  know  the  exact 
solution.  Below  we  give  some  more  physical  examples  using  material  data  for  a  Nylon  6.6 
compound. 

8.5  Physical  examples  with  “Maranyl” 

We  assume  a  synchronous  material  wherein  A  and  fi  exhibit  the  same  time  dependence 
(which  implies  a  constant  Poisson  ratio)  and  take  the  data  as  given  by  (7.13)  and  (7.14). 
For  the  single  stress  relaxation  function  we  take  E<p(t)  where, 

<?(*)  = 

i=0 

with: 


<P0 

=  0.183429971 

OL  0 

=  0.0 

<PI 

=  0.385804129 

ax 

=  53.821223820 

<P2 

=  0.430  765  899 

<*2 

=  1.592948754 

Here  the  tpi  are  dimensionless  while  the  a*  have  units  (years)-1.  These  data  are  derived 
from  experimental  creep  response  curves  for  Nylon  6.6  compound  “Maranyl”,  and  are 
taken  from  [45,  Equation  (5.39)].  (Note  that  we  have  “modernized”  the  units  by  using  the 
conversion  factor  lpsi  =  6894.8  Pa.) 

8.5.1  Example:  L-shaped  lever  arm 

We  now  consider  an  example  of  an  L-shaped  lever  arm  as  shown  in  Figure  8.2.  The  arm  is 
fixed  rigidly  in  both  displacements  along  its  leftmost  vertical  edge.  In  addition  a  vertical 
traction  of  —5  MPa  is  applied  along  the  horizontal  top  edge,  and  a  horizontal  traction  of 
—500 y  kPa  is  applied  along  the  rightmost  vertical  edge.  All  other  data  is  as  before  (e.g. 
k  =  1  and  T  =  1). 

The  first  tabulated  results  are  for  uniform  refinements  with  h  =  0.4,  0.2,  0.1, _ 


ki 

N  «l«a 

llulLE 

lluhll.E 

I  uh  1  +  i  •  1 

|uh|+««t f#| 

I lu-uhl I _E 

•«t  I  1 • 1 1 _E 

inf  I  u-uh  | 

I  I u-uh | | 

1  I «-*hl 1 

1.0000# +00 

11 

0.000#+00 

2.666#+02 

3. 629# +02 

2.690#+02 

2 . 666369#+02 

3.472313#+01 

1.784# -02 

9.893#-03 

3 . 200#+06 

1.0000# +00 

so 

0.000#+ 00 

3.080#+02 

4.356«+02 

3. 109#+02 

3,080140#+02 

4 . 257241#+01 

2.752«-02 

1 . 549#-02 

4. 162#+06 

1.0000# +00 

204 

0.000#+00 

3.643#+02 

6. 010# +02 

3. 558«+02 

3.642776#+02 

3.285800#+01 

3.701#-02 

2 . 169#-02 

4.812#+06 

1 . 0000#+00 

890 

0.000#+00 

3.693«+02 

6. 223# +02 

3.700#+02 

3.693209#+02 

2.215892«+01 

4.005#-02 

2. 329#-02 

6.017#+06 

1.0000# +00 

3744 

0.000#+00 

3.764#+02 

5.323«+02 

3 .765#+02 

3.763879#+02 

1 . 080548#+01 

4. 146#-02 

2.414#-02 

6.072#+06 

1.0000# +00 

15208 

0.000#+00 

3.781#+02 

6.347#+02 

3 . 781#+02 

3.780726#+02 

7.078166#+00 

4. 172«-02 

2. 430# -02 

6. 094# +06 
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The  next  table  of  results  are  for  adapted  solutions  with  TOL  =  32,  24,  16,  8  and 
where  we  again  start  with  an  initial  mesh  of  h  =  1.0.  Some  of  the  meshes  and  plots  of  the 
stress  surfaces  are  shown  in  the  Figure  8.2. 


ki 

N  *l*m 

llulLE 

lluhll.E 

|uhl+l*l 

|uh|+«*t |*| 

1 1 u-uh I  I _E 

•at  ||e| | _E 

inf | u-uh | 

I J u-uh I  I 

lls-Bh|| 

1.0000* +00 

308 

0.000*+00 

3. 642* +02 

5. 151*+02 

3.653*+02 

3.642289*+02 

2.733388*+01 

3.884*-02 

2. 268* -02 

4.9360+06 

1.0000* +00 

608 

0.000*+00 

3.715*+02 

6.254*+02 

3.721*+02 

3. 714864* +02 

2.074337*+01 

4.034*-02 

2.357.-02 

6.018*+06 

1.0000* +00 

1542 

0.0000+00 

3.757*+02 

5.314«+02 

3.760*+02 

3. 767406* +02 

1.441982*+01 

4. 119«-02 

2 . 402*-02 

5.0650+06 

1.0000«+00 

6460 

0.000*+00 

3.782*+02 

6.348«+02 

3.782«+02 

3.781621*+02 

7. 942375* +00 

4.170*-02 

2.430e-02 

5.0930+06 

Here  it  requires  over  15,000  elements  to  achieve  a  similar  estimated  error  as  the 


adaptive  solution  produces  with  only  6460  elements. 


8.5.2  Example:  a  simple  crack 


We  now  consider  a  simple  horizontal  crack  in  a  rectangular  component.  Due  to  symmetry 
we  consider  only  the  upper  half  of  the  component,  and  load  the  top  edge  with  a  vertical 
traction  of  5  MPa — see  Figure  8.3. 

Following  exactly  the  same  pattern  as  above,  the  first  tabulated  results  are  for  uniform 
refinements  with  h  =  0.4,  0.2,  0.1,  — 


ki 

N  *l«m 

HulLE 

lluhILE 

luhl+l* 1 

|uh|+0«t |«l 

I lu-uhl l_E 

0>t  l!*f |„E 

inf lu-uhl 

1 lu-uhl | 

1 |a-«hl 1 

1.0000* +00 

8 

0.0000+00 

2 . 130*+02 

3.0120+02 

2.1440+02 

2.129519*+02 

2.5226030+01 

7.835a-03 

2.5420-03 

2.8330+06 

1.00000+00 

34 

0.0000+00 

2 . 303*+02 

3.2560+02 

2.3200+02 

2.3026290+02 

2.8013960+01 

1.012.-02 

3.2250-03 

3.1720+06 

1.00000+00 

146 

0.0000+00 

2.4440+02 

3.4560+02 

2.4530+02 

2.444114*+02 

2.1000320+01 

1.2160-02 

3.836*-03 

3.3560+06 

1.00000+00 

618 

0.000«+00 

2.614*+02 

3.5550+02 

2.5190+02 

2.514081*+02 

1.6205630+01 

1.3090-02 

4.1390-03 

3.4600+06 

l.OOOOa+OO 

2596 

0.0000+00 

2.567*+02 

3.630*+02 

2.5690+02 

2.5670100+02 

9.7006420+00 

1.379*-02 

4.3730-03 

3.5400+06 

1.00000+00 

10394 

0.0000+00 

2.5870+02 

3.6590+02 

2.588«+02 

2.5874600+02 

6.990774«+00 

1.4050-02 

4.462*-03 

3.5730+06 

The  next  table  of  results  are  for  adapted  solutions  with  TOL  =  32,  24,  16,  8  and 
where  we  again  start  with  an  initial  mesh  of  h  =  1.0.  Some  of  the  meshes  and  plots  of  the 
stress  surfaces  are  shown  in  the  figure. 


ki 

N  *l0io 

llulLE 

lluhILE 

luhl+lo 1 

|uh|+*rt  1  •  1 

1 lu-uhl 1 _E 

•  Bt  IUILE 

inf  |  u-uh  1 

I lu-uhl | 

lli-ih|| 

1.00000+00 

146 

0.0000+00 

2.444*+02 

3.4560+02 

2.4530+02 

2.4441140+02 

2.1000320+01 

1.2160-02 

3.8360-03 

3.3560+06 

1.00000+00 

146 

0.0000+00 

2.444*+02 

3.4560+02 

2.4530+02 

2.4441140+02 

2.100032*+01 

1.2160-02 

3.8360-03 

3.3560+06 

1.00000+00 

301 

0.0000+00 

2.5430+02 

3.6970+02 

2.5480+02 

2.543484*+02 

1.5923600+01 

1.3400-02 

4.2550-03 

3.6200+06 

1.00000+00 

2235 

0.0000+00 

2.5960+02 

3.6720+02 

2.5970+02 

2.6961600+02 

7.1571020+00 

1.4130-02 

4.4960-03 

3.5910+06 

Here  it  requires  over  10,000  elements  to  achieve  a  similar  estimated  error  as  the 


adaptive  solution  produces  with  only  2235  elements. 


8.5.3  Example:  webbed  angle  bracket 

Our  next  example  is  of  a  webbed  angle  bracket  as  shown  in  Figure  8.4.  The  bracket 
is  constrained  both  horizontally  and  vertically  at  the  two  lower  horizontal  “pads”  and 
traction-loaded  only  at  the  top  vertical  pad  on  the  right.  The  horizontal  loading  is  -5  MPa 
and  the  vertical  —  500  kPa. 

We  again  begin  by  tabulating  resukts  for  uniform  refinements  with  h  =  0.4,  0.2, 

0.1,.... 
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ki  K  *1*« 

1.0000«+00  33 

1.0000«+00  39 

l.OOOOe+OO  78 
1.0000«+00  276 

1.0000*+00  1488 

1.0000*+00  6235 

1.0000e+00  25466 


llulIJS 

lluhILE 

|uhl+l«! 

|uh|+«*t|«| 

Mu-uhlLE 

•*t  IUILE 

inf  1 u-uh 1 

1 1 u-uh | | 

1 ii-ih| 1 

0.000* +00 

1 . 790*+02 

2.632«+02 

1.798*+02 

1.790362«+02 

1.699838«+01 

1.919*-02 

4.961*-03 

2 . 203*+06 

0.000* +00 

1 . 845*+02 

2. 609* +02 

1.852*+02 

1.844667* +02 

1.68 1407* +01 

2.007*-02 

6.352*-03 

2 . 305*+06 

0.000 *+00 

1.910*+02 

2.701*+02 

1.922*+02 

1.910198*+02 

2. 147025* +01 

2 . 115* -02 

6.009*-03 

2. 479* +06 

0.000*+00 

2.082e+02 

2.944.+02 

2.091e+02 

2. 081676* +02 

1.986391*+01 

2. 438* -02 

7.869*-03 

2.717*+06 

0.000* +00 

2.1906+02 

3.097*+02 

2 . 193*+02 

2.190010.+02 

1.118244*+01 

2. 800* -02 

9.279*-03 

2.869«+06 

0.000*+00 

2. 220* +02 

3. 139*+02 

2.221«+02 

2.219738«+02 

6.972105*+00 

2. 910* -02 

9.670*-03 

2. 913* +06 

0.000«+00 

2.231*+02 

3. 155*+02 

2.2316+02 

2.2307826+02 

4.052005*+00 

2. 948* -02 

9.802*-03 

2.9256+06 

The  next  table  of  results  are  for  adapted  solutions  with  TOL  =  20,  15,  10,  5  and  where 
we  now  start  with  an  initial  mesh  of  h  —  0.05  (due  to  the  slender  nature  of  the  component). 
Some  of  the  meshes  and  plots  of  the  stress  surfaces  are  shown  in  the  Figure  8.4. 


ki  K  «1mi  llull.E  lluhILE  |uh|+l*l  |uh|+«it|*|  |Ju-ohll_E  «it  1MI.E  inf|u-uh|  | |u-uhl I  ||i-sh|| 

1.0000*+00  276  0.000«+00  2.082*+02  2.944«+02  2.091«+02  2.081676*+02  1.986391*+01  2.438«-02  7.869«-03  2.717«+06 

1.0000«+00  829  0.000*+00  2.179«+02  3.082«+02  2.183«+02  2.179416*+02  1.306965«+01  2.768«-02  9.049«-03  2.859«+06 

1.0000«+00  2487  0.000«+00  2.213«+02  3.130«+02  2.215*+02  2.213019e+02  9.014411«+00  2.882«-02  9.529«-03  2.903«+06 

1.0000«+00  12232  0,000«+00  2.230a+02  3.154*+02  2.231«+02  2.230293«+02  4.942654«+00  2.945a-02  9.777*-03  2.925«+06 


Here  it  requires  over  25,000  elements  to  achieve  a  similar  estimated  error  as  the 
adaptive  solution  produces  with  only  12,232  elements. 

Each  of  these  examples  clearly  demonstrates  how  adaptive  mesh  refinement,  guided 
by  an  a  posteriori  error  estimate,  can  result  in  an  acceptable  solution  requiring  far  fewer 
elements  that  would  be  needed  when  using  uniform  meshes. 


sigmall 
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Figure  8.1:  Adapted  meshes  for  the  exact  solution  with  TOL  =  15,  10  and  5  (initial  mesh: 
h  =  0.1).  Also  shown  are  plots  of  the  stress  surfaces  for  TOL  =  10. 


391  nodes,  695  elements 
85  boundary  sides  during 
times  t  G  (0.000000, 1.000000). 


595  nodes,  1056  elements 
132  boundary  sides  during 
times  t  e  (0.000000, 1.000000) 


1656  nodes,  3064  elements 
246  boundary  sides  during 
times  t  G  (0.000000,1.000000). 


sigmal  1  at  time  =  1 .000000 


sigma12  at  time  =  1.000000 


sigma22  at  time  =  1.000000 
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Figure  8.2:  Adapted  meshes  for  the  L-shaped  lever  arm  with  TOL  =  32,  24  and  16  (initial 
mesh:  h  =  0.1).  Also  shown  are  plots  of  the  stress  surfaces  for  TOL  =  8. 


186  nodes,  308  elements 
62  boundary  sides  during 
times  t  £  (0.000000, 1.000000). 


347  nodes,  608  elements 
84  boundary  sides  dining 
times  t  £  (0.000000,1.000000). 


836  nodes,  1542  elements 
128  boundary  sides  during 
times  t  £  (0.000000,1.000000). 


sigmall  at  time  =  1.000000 


sigma12  at  time  =  1 .000000 


0  o 


sigma22  at  time  =  1 .000000 
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Figure  8.3:  Adapted  meshes  for  the  simple  symmetric  crack  with  TOL  =  24,  32  and  8 
(initial  mesh:  h  =  0.1).  Also  shown  are  plots  of  the  stress  surfaces  for  TOL  =  8. 


89  nodes,  146  elements 
30  boundary  sides  during 
times  t  £  (0.000000, 1.000000). 


1191  nodes,  2235  elements 
145  boundary  sides  during 
times  t  £  (0.000000, 1.000000). 


sigmal  1  at  time  =  1.000000 


174  nodes,  301  elements 
45  boundary  sides  during 
times  t  £  (0.000000, 1.000000). 


sigmal  2  at  time  =  1.000000 


0  o 
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Figure  8.4:  Adapted  meshes  for  the  webbed  angle  bracket  with  TOL  =  15,  10  and  5  (initial 
mesh:  h  =  0.05).  Also  shown  are  plots  of  the  stress  surfaces  for  TOL  =  10. 


507  nodes,  829  elements 
185  boundary  sides  during 
times  t  E  (0.000000, 1.000000). 


1393  nodes,  2487  elements 
299  boundary  sides  during 
times  t  €  (0.000000, 1.000000). 


6445  nodes,  12232  elements 
658  boundary  sides  during 
times  t  E  (0.000000, 1.000000) 


sigmall  at  time  =  1.000000 


sigma  1 2  at  time  =  1 .000000 


sigma22  at  time  =  1.000000 
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Part  III 
Closure 
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Chapter  9 


Obtaining  and  using  the  software 


9.1  Introduction 

The  C  code  used  to  generate  the  numerical  solutions  detailed  in  the  previous  chapters 
is  available  via  ftp.  This  chapter  describes  how  to  retrieve  the  source  files  and  compile 
and  run  them  on  an  X- window  Unix  platform.  We  also  give  a  short  “manual”  describing 
how  the  code  can  be  set  up  to  solve  a  specific  problem.  We  do  this  by  choosing  a  simple 
example  problem  and  demonstrating  how  the  input  files  should  be  generated  to  model  the 
domain  and  boundary  conditions.  Since  the  software  comes  packaged  with  the  example 
problems  used  earlier  in  this  report,  these  together  with  this  manual  should  provide  all 
the  necessary  information. 

Note:  this  software  is  what  we  term  a  “research  code”.  It  has  been  developed  in 
a  piecemeal  fashion  over  a  long  period  of  time  as  and  when  new  research  results  become 
available  and  are  suitable  for  implementation.  Consequently,  no  claims  are  made  for  it 
being  efficient  or  robust,  and  it  should  not  be  used  in  a  situation  where  its  output  may 
have  a  safety  or  financial  implication. 


9.2  Obtaining  the  software 


Firstly,  connect  to  the  Brunei  University  ftp  server  via  the  command, 
ftp  ftp.brunel.ac.uk 


When  prompted  logon  as  “anonymous”  and  enter  your  email  address  as  a  password.  When 
you  get  the  prompt  issue  the  following  commands: 


cd  icsrsss 
bin 

get  seedproj.tar 
quit 
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You  should  now  have  a  binary  file  called  seedproj  .tar  in  your  local  working  directory. 
This  is  a  tape-archived  record  of  the  source  codes  created  using  the  Unix  tar  utility. 

The  next  step  is  to  extract  the  directory  structure.  Do  this  with  tar  using  the 
command, 

tar  xvf  seedproj. tar 

You  now  should  have  a  directory  called  seedproj  which  contains  the  sub-directories  hold¬ 
ing  the  source  codes  for  the  mesh  generator  and  finite  element  solver.  At  this  point  the 
archive  file  seedproj  .tar  can  be  deleted. 

To  compile  the  mesh  generator  issue  the  following  commands: 

cd  seedproj /advance/source 
make. 

and  to  compile  the  solver  type, 

cd  . . / . . /quasivis/native/source 
make 

There  is  also  a  very  simple  X- window  plotting  program  called  Xwin  which  will  be  useful 
later  when  generating  the  mesh.  To  compile  this  type, 

cd  /Xwin 

make 

If  this  stage  has  gone  well  you  should  have  three  executable  files,  advance,  fern  and  Xwin 
in  the  directory, 

seedproj /bin 

Now  change  to  the  directory  seedproj  /quasi  vis /data  and  list  the  contents.  You 
will  see  subdirectories  such  as  leverarm,  crack,  etc.  These  are  the  input  files  for  the 
numerical  examples  shown  in  previous  chapters.  The  source  files  are  supplied  and  set  up 
so  as  to  solve  the  problem  with  an  exact  solution  (see  previously  in  Subsection  8.4.2)  and 
to  check  that  the  installation  was  successful  type, 

adapt  exact  15 

This  executes  the  sh  script  file  adapt  and  solves  the  problem,  with  TOL  =  15,  as  specified 
by  the  files  in  the  exact  directory.  Specifically  the  mesh  generator  advance  builds  the 
mesh  according  to  the  files  in  the  exact  directory,  and  then  adapt  launches  the  finite 
element  solver  fern.  After  the  make  (i.e.  compilation)  stage  you  should  see  an  X-window 
popped  to  the  screen  showing  a  sequence  of  adapted  meshes.  You  will  also  see  many  lines 
of  text  rolling  around  the  terminal  screen.  The  end  of  the  textual  output  should  look  like 
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ki  H  «l«n  I  lul  |_E  lluhll.E  |uh|  +  M  |uh|+e«t|e|  |  |u-uh|  |_E  «st  tl«li_E  influ-uhl  ilu-uh||  ||«-«h|f 

1.0000e+00  695  8.398e+01  8.331e+01  8.398e+01  8.416e+01  1.061322e+01  1.186403e+01  2.016e-05  3.930e-06  4.303e+0S 

Tu«  Dec  1  17:15:34  GMT  1998 

(although  you  will  certainly  have  a  different  time  stamp  showing).  Note  that  these  figures 
agree  with  those  shown  in  the  first  line  of  the  table  of  adapted  solutions  in  Subsection  8.4.2. 
If  all  of  this  happens  the  installation  went  well,  and  you  may  now  define  your  own  domain 
and  quasistatic  viscoelasticity  problem. 


9.3  Defining  the  domain 

As  the  name  advance  suggests,  the  mesh  generator  is  a  (simple)  implementation  of  the 
advancing  front  technique.  For  this  it  is  necessary  only  to  define  the  boundary  of  the 
domain  according  to  a  simple  orientation  rule,  and  provide  a  function  h  =  h(x ,  y)  that 
describes  the  desired  mesh  size  variation  over  the  domain.  In  our  implementation  h  = 
h(x,y)  is  specified  as  a  constant  ho:  which  we  term  the  basic  mesh  size,  and  local  mesh 
features  and  size  control  are  effected  through  source  points  (in  the  file  srcpnts.dat), 
source  lines  (in  the  file  srcline.dat),  source  circles  (in  the  file  srccirc.dat)  and 
source  discs  (in  the  file  srcdisc.dat).  We’ll  get  to  these  later  but  for  now  note  that  in 
this  way  an  a  priori  graded  mesh  can  be  created  by  suitably  editing  one  of  these  input 
files  rather  than  hard  coding  a  C  function,  double  h  (double  x,  double  y),  and  then 
having  to  re-compile  each  time. 

At  the  moment  only  two-dimensional  polygonal  domains  are  supported  by  advance, 
with  the  boundary  defined  piecewise  by  the  endpoints  of  straight-line  segments.  The 
domain  can  however  be  multiply-connected  (i.e.  have  an  arbitrary  number  of  “holes”), 
and  it  would  not  be  difficult  in  the  future  to  incorporate  curved  boundaries. 

To  illustrate  how  to  set  up  a  domain  in  the  way  required  for  advance  we  will  consider 
the  case  of  a  trestle-shaped  structure  with — for  illustration  purposes  only — a  hole  in  it. 
The  domain  is  shown  in  Figure  9.1  and  is  supposed  to  symmetrical  about  the  vertical  edge 
at  x  =  11cm.  The  lowest  horizontal  edge  is  constrained  to  have  no  vertical  displacement 
and,  to  impose  the  symmetry,  the  right-most  vertical  edge  is  constrained  to  have  no 
horizontal  displacement.  We  assume  that  a  load  (i.e.  surface  traction)  of  <72  —  —  1  MPa 
acts  vertically  downward  on  the  uppermost  horizontal  edge. 

Each  corner  of  the  trestle  has  a  numbered  node  associated  with  it  and  in  between 
each  node  we  identify  a  numbered  (shown  in  bold)  boundary  edge.  Note  that  the  node  and 
edge  numbering  must  begin  at  zero  (this  is  because  all  arrays  begin  with  a  zero  subscript 
in  the  C  language).  Also,  although  the  node  and  edge  numbering  follow  a  simple  counter¬ 
clockwise  pattern  in  this  example,  there  is  absolutely  no  need  for  this.  The  crucial  thing 
is  that  the  orientation  rule  which  we  describe  below  is  observed  when  defining  the  edge 
connectivity. 

Before  constructing  the  input  files  we  need  to  assign  the  correct  boundary  condition 
codes  to  the  nodes  and  edges.  These  integer  values  describe  the  type  of  constraint  (if 
any)  that  the  node  or  edge  is  subject  to,  and  are  given  by: 


0  —  unfixed  (for  interior  nodes  only). 
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Figure  9.1:  The  left  half  of  the  trestle  domain.  (Measurements  are  in  cm.) 


8  9  10  11  12 


1  —  Non-homogeneous  Dirichlet  condition  -  unused  at  the  moment. 

2  —  Homogeneous  Dirichlet  condition  (zero  displacement). 

3  —  Non-homogeneous  Neumann  condition  (prescribed  traction). 

4  —  Homogeneous  Neumann  condition  (zero  traction). 

Each  node  and  each  edge  is  now  assigned  a  triple  of  integers  made  up  from  the  above  codes 
(except  for  0  —  unfixed,  this  is  reserved  for  interior  nodes  created  by  the  mesh  generator). 
The  first  two  elements  of  the  triple  correspond  to  displacements  or  tractions  in  the  x  and 
y  direction  respectively,  while  the  third  element  is  unused  at  present  (and  so  its  value  is 
irrelevant). 

This  works  in  the  following  way.  Nodes  3  and  4,  and  edge  3,  at  the  foot  of  the  trestle 
are  to  be  constrained  so  as  to  have  no  vertical  displacement  and  this  is  accomplished  by 
setting  a  2  in  the  second  element  of  the  triple.  Also,  we’ll  assume  that  there  is  no  traction 
acting  horizontally  along  edge  3  and  so  the  first  element  of  the  triple  is  set  as  4.  For  this 
part  of  the  boundary  we  then  get, 

boundary  codes  for  node  3  are:  4  2  0 
boundary  codes  for  node  4  are:  4  2  0 
boundary  codes  for  edge  3  are:  4  2  0 

Now,  for  edge  0,  defined  by  nodes  0  and  1,  we  assume  zero  horizontal  traction  and  this 
leads  to  us  setting  the  first  element  of  the  triple  as  4.  On  the  other  hand  we  have  a 
non-zero  vertical  traction  acting  on  this  edge  and  so  the  second  element  is  set  as  3.  This 
gives: 
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boundary  codes  for  edge  0  are:  4  3  0 

By  symmetry  node  0  is  fixed  horizontally  -  code  2  -  but  is  subject  to  vertical  traction  - 
code  3.  Thus: 

boundary  codes  for  node  0  are:  2  3  0 
boundary  codes  for  node  1  are:  4  3  0 

since  node  1  is  free  and  unforced  in  the  x  direction. 

Now  let  us  take  a  look  at  edge  6  as  defined  by  nodes  6  and  0.  Due  to  symmetry  this 
edge  may  not  move  horizontally  -  code  2  -  and  must  be  free  of  vertical  traction  -  code  4. 
We  get  the  same  for  node  6  and  so: 

boundary  codes  for  node  6  are:  2  4  0 
boundary  codes  for  edge  6  are:  2  4  0 

Every  other  node  and  edge  is  free  to  move  (so  we  cannot  use  either  of  codes  1  or  2),  but  is 
also  traction-free.  Thus  for  all  other  nodes  and  edges  we  set  the  boundary  code  as:  4  4  0. 

Equipped  with  this  information  we  now  create  the  first  input  files  for  the  mesh  gen¬ 
erator.  Inside  the  seedproj/quasivis/data  directory  create  a  new  subdirectory  called 
trestle.  Change  into  this  new  directory  and  edit  the  files  def start .  inf,  def nodes .  inf 
and  def  sides,  inf  as  shown  below  (do  not  include  the  underlined  headings). 

The  input  file  c def start .inf  * 


Number _of_nodes:  11 
Number_of_sides:  11 

bounding_box_f orographies  0.0  0.12  -0.03  0.09 

The  file  def  start .  inf  simply  shows  the  number  of  nodes  and  (necessarily  the  same 
for  polygonal  domains)  edges  used  to  define  the  domain — in  this  case  11.  The  file  also 
defines  a  bounding  box  to  control  the  graphical  display  of  the  domain.  This  box  is 
defined  by  a  line  of  the  form, 

bounding_box  jf  or  _gr  aphi  c  s  x\ow  ^high  2/low  2/high 

which  describes  the  box  by  the  diagonal  running  from  (xiOWJ2/iow)  to  (^high, 2/high)*  The 
bounding  box  is  assumed  to  be  square  by  all  of  the  graphics  routines  and  so  it  is  necessary 
to  ensure  that 


S'high  #low  —  2/high  2/low  • 


If  this  does  not  hold  then  the  display  will  be  distorted  by  a  relative  scale  difference  in  the 
horizontal  and  vertical  directions.  For  our  example  we  will  specify  all  lengths  in  metres 
and  so  we  choose  a  bounding  box  with  the  line, 
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bounding_box_f orographies  0.0  0.12  -0.03  0.09 

The  domain  fits  inside  this  bounding  box  and  so  should  appear  centered  in  the  graphical 
displays.  By  altering  the  bounding  box  one  can  scale  the  x  and  y  dimensions  of  the  domain 
in  different  ways  or  place  it  off-center,  as  well  as  reducing  its  apparent  distance  from  the 
eye. 

The  input  file  ‘defnodes . inf 3 


0 

2 

3  0 

0.11 

0.06 

1 

4 

3  0 

0.01 

0.06 

2 

4 

4  0 

0.01 

0.04 

3 

4 

2  0 

0.02 

0.01 

4 

4 

2  0 

0.04 

0.01 

5 

4 

4  0 

0.05 

0.04 

6 

2 

4  0 

0.11 

0.04 

7 

4 

4  0 

0.04 

0.04 

8 

4 

4  0 

0.02 

0.04 

9 

4 

4  0 

0.02 

0.05 

10 

4 

4  0 

0.04 

0.05 

The  file  def nodes. inf  contains  the  numbered  (in  order)  list  of  nodes  defining  the 
domain,  with  each  integer  label  followed  by  the  boundary  code  triple  and  the  node  coordi¬ 
nates.  There  is  nothing  special  about  the  formatting  or  spacing  of  the  fields  on  each  line, 
so  long  as  at  least  one  whitespace  character  separates  each. 

The  input  file  ‘def sides . inf ’ 


0  4  3  0  0  1 

1  4  4  0  1  2 

2  4  4  0  2  3 

3  4  2  0  3  4 

4  4  4  0  4  5 

5  4  4  0  5  6 

6  2  4  0  6  0 

7  4  4  0  7  8 

8  4  4  0  8  9 

9440  9  10 

10  4  4  0  10  7 


The  file  def  sides,  inf  is  similar  to  defnodes.inf;  it  contains  the  numbered  list  of 
edges  defining  the  domain.  On  each  line  the  integer  edge  label  is  followed  by  the  boundary 
condition  triple,  and  then  by  the  pair  of  nodes  that  define  the  edge.  For  example,  the  line 


3  4  2  0  3  4 
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shows  that  edge  3  takes  the  boundary  codes  4  2  0  and  is  defined  by  nodes  3  and  4.  The 
order  in  which  these  two  nodes  appear  is  crucial  and  must  obey  the  orientation  rule. 
This  rule  is  simple  enough:  it  says  that  in  travelling  along  the  edge  from  the  first  node  to 
the  second  node  the  domain  must  be  on  the  left.  For  this  reason  all  the  edges  on  the  outer 
boundary  (i.e.  0  to  6)  are  described  by  a  counter-clockwise  node  ordering,  while  all  the 
edges  on  the  interior  boundary  (i.e.  edges  7  to  10)  are  described  by  a  clockwise  ordering. 


9.4  Grading  and  generating  the  mesh 


As  far  as  the  mesh  generator  advance  is  concerned  the  domain  is  now  completely  defined, 
but  we  cannot  yet  use  it  to  create  the  mesh  because  we  have  to  specify  the  mesh  size 
function  h  =  h(x,y).  This  function  should  return  the  desired  size  of  the  triangles  in 
the  vicinity  of  the  point  (x7y).  To  do  this  we  need  the  files  srcpnts.dat,  srcline.dat, 
srccirc.dat  and  srcdisc.dat.  These  files  can  be  used  to  specify  the  basic  value  of  h 
throughout  the  domain,  as  well  as  any  local  refinements  that  are  a  priori  required.  You 
may  find  it  easier  to  simply  copy  these  files  over  from  another  directory  (e.g.  . .  /crack) 
since  they  will  contain  helpful  annotations  at  the  bottom  that  describe  the  meaning  of  the 
information  contained  in  them. 

Once  you  have  copied  them  across  edit  the  files  srcline.dat,  srccirc.dat  and 
srcdisc.dat  and  make  sure  that  the  very  first  entry  on  the  first  line  is  the  integer  0. 
We’ll  return  to  these  files  later  but  for  now  concentrate  on  the  most  important  one: 
srcpnts.dat. 

The  first  entry  in  srcpnts.dat  is  ho?  the  basic  value  of  h  to  be  used  throughout  the 
domain.  Looking  back  at  Figure  9.1  it  seems  reasonable  to  require  an  initial  mesh  consist¬ 
ing  of  triangles  with  side  length  ho  =  0.005  metres.  To  obtain  this  we  edit  srcpnts.dat 
so  that  the  first  line  contains  the  value  0.005,  and  ensure  that  the  second  line  contains 
the  integer  0.  We  may  now  generate  the  mesh.  To  do  this  issue  the  command  (in  the 
directory  trestle), 


/../.. /bin/advance  |  /bin/Xwin 


This  initiates  the  mesh  generator  advance  and  pipes  the  output  into  the  simple  X-window 
utility  Xwin.  In  this  case  an  X-window  should  pop  to  the  screen  and  the  element  edges 
are  drawn  in  this  window  as  they  are  created.  When  the  mesh  is  complete  place  the 
mouse  cursor  in  the  X-window  and  type  the  q  button.  The  programs  will  quit  and  the 
prompt  will  return  in  the  terminal  window.  (Note  that  Xwin  is  a  very  primitive  X-window 
application  and  so  the  terminal  may  display  “junk”  output  while  it  is  running.) 

If  you  now  list  the  files  in  the  trestle  subdirectory  you  should  see  something  like: 


bdystart .fem 
boundary .fem 
conf ig. inf 
def nodes . inf 


def sides. inf 
def  st art  .inf 
element .fem 
memory . dat 


meshinfo.fem 
mon . out 
neighbor .fem 
nodes .fem 


srccirc.dat 

srcdisc.dat 

srcline.dat 

srcpnts.dat 
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All  except  the  def  * .  inf  and  src*  .dat  files  that  we  created  earlier  are  disposable  (in  that 
they  can  be  regenerated),  and  so  can  eventually  be  removed  with  a  command  line  like, 

rm  * . fern  memory.dat  mon.out  scratch.dat 

(Note  that  the  files  mon.out  and  scratch.dat  may  or  may  not  be  present,  depending  on 
how  advance  terminated.  These  files  are  always  disposable  and  may  use  up  a  lot  of  disk 
space — it  is  advisable  therefore  to  always  remove  them.)  The  files  *.fem  are  used  by  the 
finite  element  code  to  determine  the  initial  mesh.  Now  we  have  these  *.fem  files  we  can 
set  up  and  execute  the  finite  element  calculation  in  the  next  section,  but  for  the  remainder 
of  this  section  we’ll  explain  a  little  more  about  how  to  use  the  src*  .dat  files  to  control  a 
priori  mesh  grading. 

The  mesh  that  appeared  on  the  screen  should  have  looked  like  that  on  the  left  in 
Figure  9.2  and,  one  might  think,  is  perfectly  acceptable  for  the  initial  mesh  in  an  adaptive 
calculation.  However,  it  is  also  true  that  an  adaptive  calculation  can  be  much  improved  by 
a  sensible  choice  of  initial  mesh.  For  example,  for  the  trestle  it  is  likely  that  high  stresses 
will  be  present  at  the  re-entrant  corner  at  (0.05,0.04),  and  so  it  is  sensible  to  build  this 
information  into  the  initial  mesh  by  asking  for  a  finer  mesh  grading  at  this  point.  This 
can  done  by  introducing  a  source  point  at  this  corner  by  making  an  entry  in  the  file 
srcpnts.dat.  Edit  this  file  so  that  the  first  three  lines  are  as  shown  below. 

0.005 

1 

0.05  0.04  0.02  0.002 

Recall  that  the  first  entry  is  no  more  than  ho  =  0.005m.  The  second  entry  states 
that  the  mesh  is  to  contain  one  source  point,  and  the  third  line  then  defines  this  source 
point  to  be  located  at  (0.05,0.04),  to  have  a  radius  of  influence  of  0.02m  and  to  have  a 
local  mesh  size  at  (0.05,0.04)  of  0.002m.  The  mesh  generator  will  now  attempt  to  create  a 
mesh  with  this  local  size  at  the  re-entrant  corner,  blending  smoothly  to  the  basic  mesh  size 
throughout  a  circle  of  radius  0.02m.  Run  the  mesh  generator  again  with  the  command, 

../../. ./bin/advance  I  /bin/Xwin 

and  you  should  now  get  the  mesh  shown  on  the  right  of  Figure  9.2. 

On  an  intuitive  level  this  seems  a  much  more  sensible  choice  for  the  initial  mesh,  and 
it  is  simple  to  specify  through  the  use  of  the  input  file  srcpnts.dat.  Any  number  of 
source  points  can  be  added  in  this  way,  with  the  total  number  given  in  the  second  line  of 
the  file,  so  long  as  two  simple  rules  are  followed.  The  local  mesh  size  should  always  be 
smaller  than  the  basic  mesh  size,  and  the  radius  of  influence  should  be  at  least  twice  the 
basic  mesh  size. 

Source  lines,  source  circles  and  source  discs  can  be  specified  in  a  similar  way  through 
the  input  files  srcline.dat,  srccirc.dat  and  srcdisc.dat.  To  illustrate  this  edit  the 
file  srcline .  dat  so  that  the  first  five  lines  are  given  by, 
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Figure  9.2:  Trestle  mesh  generated  by  advance.  The  mesh  on  the  left  has  a  basic  value  of 
ho  =  0.005m  throughout  while  the  one  on  the  right  has  a  source  point  at  (0.05, 0.04). 


a 

H 


0.04  0.04  0.02  0.04 
0.02  0.04  0.02  0.05 
0.02  0.05  0.04  0.05 
0.04  0.05  0.04  0.04 


0.005  0.002 
0.015  0.002 
0.015  0.002 
0.005  0.002 


This  now  specifies  a  source  line  for  each  edge  of  the  rectangular  “hole”  in  the  trestle.  The 
annotations  in  the  file  explain  what  the  entries  on  each  line  actually  mean,  and  running 
the  mesh  generator  again  gives  the  mesh  shown  on  the  left  of  Figure  9.3. 

The  mesh  now  appears  to  be  suitable  for  input  to  the  finite  element  code  and  indeed 
this  is  the  one  we  will  use  below.  However,  we  will  give  one  more  example  of  how  to  achieve 
local  mesh  grading.  To  create  a  source  circle  of  radius  0.008m  at  the  point  (0.08,0.05), 
with  a  local  mesh  size  0.001m  and  radius  of  influence  0.01m,  edit  the  file  srccirc.dat  so 
that  the  first  two  lines  are: 

1 

0.08  0.05  0.008  0.01  0.001 

The  first  line  specifies  how  many  source  circles  there  are  while  the  subsequent  lines  (in 
this  only  one)  specify  the  circles  themselves.  The  resulting  mesh  is  shown  on  the  right  of 
Figure  9.3.  We  have  included  this  source  circle  for  demonstration  purposes  only,  it  will  be 
of  no  further  use  so  change  the  1  in  the  first  line  in  the  file  back  to  0  to  switch  the  feature 
off. 


9.5  The  finite  element  calculation 


We  now  have  our  initial  mesh,  the  one  on  the  left  of  Figure  9.3,  and  so  we  are  ready  to  feed 
it  in  to  the  finite  element  code.  However,  to  drive  the  solver  we  must  first  provide  the  file 
config.inf  in  the  trestle  directory.  This  file  contains  basic  configuration  parameters 
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Figure  9.3:  Trestle  mesh  generated  by  advance.  The  mesh  on  the  left  has  source  lines 
around  the  perimeter  of  the  square  hole  while  the  one  on  the  right  shows  the  effect  of  a 
source  circle. 


for  the  code  and,  again,  is  best  created  by  first  copying  an  existing  version  over  from 
another  directory  (e.g.  . .  /crack)  and  then  editing  it  so  that  it  appears  as  below. 

1.0  =  k_init  (used  only  if  not  set  on  command  line) 

1.0  =  T,  the  final  time  in  the  time  interval  (0,T) 

1.0e-7  Conjugate  gradient  residual  tolerance 

1.0e-2  Error  tolerance,  non-specific  at  the  moment 

5  5  400  400  graphics  window:  x,  y,  width,  height 

2.0  magnification  scale  for  graphical  output 

15  (or  25)  minimum  permitted  angle  in  the  mesh  (in  degrees) 

150  (or  120)  maximum  permitted  angle  in  the  mesh  (in  degrees) 

5  "memscale"  controls  memory  allocation  for  the  mesh 


This  file  contains  algorithm  parameters. 

Note  that  only  the  numeric  data  appearing  on  the  left  of  each  line  are  used,  the  text  to 
the  right  being  a  short  description  of  what  the  data  are.  We’ll  describe  the  meaning  of 
each  of  these  in  a  little  more  detail. 

k_init:  This  is  the  initial  time  step — it  can  be  overridden  by  a  command  line  option. 

T:  The  final  time  to  calculate  to  where  the  time  interval  is  assumed  finite:  J  :=  [0,T]. 

Since  the  code  is  not  yet  ready  to  perform  adaptive  time  stepping  we’ll  assume  a 
calculation  over  the  short  time  interval  [0,1]  in  a  single  time  step  fcjnjt  =  T  =  1.  The 
remaining  items  in  this  file  are  described  below. 

1.0e-7:  The  residual  tolerance  to  use  as  a  stopping  condition  in  the  conjugate  gradient 
solver. 
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1 .  Oe-2:  This  field  is  unused 

5  5  400  400:  These  integers  control  the  location  and  size  of  the  X-window  that  pops 
onto  the  screen.  In  this  example  the  window  will  have  its  top  and  left  hand  edges 
each  5  pixels  from  the  edge  of  the  screen,  and  will  be  400  pixels  square. 

2.0:  In  the  graphical  display  of  results  the  displacements  are  shown  on  the  mesh  magnified 
by  this  factor. 

15:  The  minimum  tolerable  interior  angle  in  a  triangle — used  by  the  mesh  adapting  rou¬ 
tines  to  determine  whether  a  triangle  should  undergo  a  “red”  or  “green”  refinement. 

150:  The  maximum  tolerable  interior  angle  in  a  triangle — as  above. 

5  "memscale":  When  the  solver  begins  to  execute  it  needs  to  set  aside  enough  memory 
to  store  the  mesh.  However,  in  an  adaptive  calculation  it  is  not  known  a  priori  how 
much  memory  the  final  mesh  will  require.  Thus  a  mesh  memscale  times  denser  than 
the  initial  mesh  is  assumed.  This  field  can  be  overridden  on  the  command  line  as  we 
show  below,  and  in  future  versions  of  the  code  we  will  probably  implement  a  mesh 
reallocation  routine  so  as  to  exploit  the  ability  of  C  to  perform  dynamic  memory 
allocation.  This  parameter  will  then  be  redundant. 


Now  we  come  to  the  least  user-friendly  part  of  the  set-up,  specifying  the  loads.  Change 
into  the  directory  seedpro j  /quasivis/native/source.  Here  you  will  find  all  the  C  source 
code  for  the  finite  element  solver.  The  file  that  we  need  to  edit  to  specify  the  loads  is 
lesol .  c.  Unfortunately  a  moderate  amount  of  C  coding  is  required  here,  and  below  we’ll 
go  over  the  main  points.  Load  lesol. c  into  an  editor  and  look  at  the  top  of  the  file. . . 

#include  "always.h" 

#include  "femdata.h" 

extern  Prony  lambda,  mu; 

/*/ 

/*  Test  exact  solutions  for  a  linear  elasticity  problem 
/*  Possible  solutions  are  selected  by  an  appropriate  #define. 

/*  The  options  are: 

/* 

/*  #define  EXACT 
/*  #define  LEVERARM 
/*  #def ine  CRACK 
/*  #def ine  WEB 
/*  #define  TRESTLE 
/*/ 

#def ine  EXACT  /*  . . .  MAKE  THE  CHOICE  */ 


trigonometric  displacements  -  see  notes, 
actual  problem,  nothing  known 
actual  problem,  nothing  known 
actual  problem,  nothing  known 
actual  problem,  nothing  known 
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The  first  #include  statements  are  directives  to  the  compiler  to  include  the  contents  of 
the  header  files  always.h  and  femdata.h.  These  can  also  be  found  in  the  directory 
native/source  and  themselves  contain  other  #include  directives.  They  need  not  worry 
us.  The  extern  statement  that  follows  next  is  there  merely  to  tell  the  compiler  that  two 
variables,  lambda  and  mu,  of  type  Prony  are  referenced  in  this  file  but  are  actually  defined 
elsewhere  (in  the  file  viscous.c  in  fact).  This  need  not  trouble  us  either.  The  next  few 
lines,  those  beginning  with  /*  are  comments  and  are  ignored  by  the  compiler.  This  is 
because  any  material  enclosed  between  an  opening  /*  and  a  closing  */  is  a  comment  in  C 
and  therefore  plays  no  role  in  the  executable  code.  Thus,  the  next  non-trivial  line  in  the 
file,  and  the  one  of  concern  to  us,  is 

#def ine  EXACT  /*  . . .  MAKE  THE  CHOICE  */ 

This  line  tells  the  compiler  that  the  string  EXACT  has  a  special  significance  in  the  source 
code  that,  follows.  As  the  comments  above  this  line  suggest,  and  with  regard  to  the 
numerical  results  presented  earlier,  it  is  also  possible  to  substitute  EXACT  with: 

EXACT  for  the  exact  solution; 

LEVERARM  for  the  L-shaped  domain; 

CRACK  for  the  simple  crack  problem; 

WEB  for  the  webbed  bracket. 

TRESTLE  for  our  trestle  problem. 

In  fact,  once  the  role  this  #def  ined  string  plays  in  the  following  code  is  appreciated,  you 
may  #def  ine  any  string  of  your  choice  in  order  to  specify  your  problem.  To  solve  our 
example  trestle  problem  change  this  line  to: 

#def ine  TRESTLE  /*  . . .  MAKE  THE  CHOICE  */ 

The  next  part  of  the  file  lesol.  c  can  be  ignored  now  until  you  come  to  the  lines 


/*** ************************************************************ ********/ 
/*  Displacements  */ 

/***********************************************************************/ 
/*  horizontal  displacements  */ 

real  uxyt(x,y,t) 
real  x,y,t; 

{ 

real  total  =  ZERO; 

#if  defined (EXACT) 

total  =  X(x,y)  *  TIME(t); 

#elif  defined (LEVERARM) 
total  =  (real)O.O; 

#elif  defined (CRACK) 
total  =  (real)O.O; 
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#elif  defined(WEB) 

total  =  (real) 0.0; 

#elif  def ined (TRESTLE) 
total  =  (real)O.O; 

#else 

fprintf (stderr,H\n\t  No  solution  def ined\n\n") ; 
#endif 

return  total; 

> 


This  declares  uxyt  to  be  a  C  function  that  takes  three  real  numbers  rr,  y  and  t,  giving 
position  and  time,  and  returns,  via  the  return  statement,  a  real  number — held  in  the  real 
variable  total.  (C  programmers  note:  the  type  real  is  used  here,  via  typedef ,  as  a  syn¬ 
onym  for  double.)  This  function  uxyt  should  return  the  horizontal  displacement  u(x ,  y,  t) 
when  it  is  known  (for  the  EXACT  solution  only)  or  zero  (for  all  the  others:  LEVERARM, 
CRACK,  WEB,  TRESTLE).  Notice  how  the  #def  ined  token  EXACT  is  used  by  the  compiler  here 
in  conjunction  with  the  “if  defined”  (#if  def  ined)  and  “else  if  defined”  (#elif  defined) 
directives.  The  meaning  and  use  of  these  is  intuitive:  since  we  had  already  #def  ined  the 
token  EXACT  then  the  compiler  only  built  the  line 

total  =  X(x,y)  *  TIME(t) ; 

into  the  executable.  (The  symbols  X  and  TIME  are  themselves  defined  further  up  the  code 
and  need  not  worry  us).  The  variable  total  is  then  set  to  a  real  value  and  returned  as 
such  every  time  the  function  is  called.  The  important  point  to  realize  here  is  that  the 
choice  of  which  branch  of  the  #if  —  #elif  to  build  in  is  made  at  compile  time,  and 
so  if  the  #defined  token  is  changed  the  code  must  be  re-compiled  to  reflect  this.  Thus 
our  newly  #defined  token  TRESTLE  requires  that  lesol.c  be  re-compiled  and  re-linked 
to  produce  an  updated  executable — this  will  happen  automatically  when  we  execute  the 
adapt  script  below. 

As  we  look  further  down  the  source  code  we  will  see  many  other  functions,  each 
constructed  in  a  similar  way.  The  “listing”  of  the  function  uxytO  given  above  shows 
how  we  include  a  new  TRESTLE  branch  into  the  #if  —  #elif  directive  and  this  must  be 
repeated  for  every  other  function  defined  the  following  list. 

uxyt  returns  real  value  of  horizontal  displacement  u(rr,y,£)  or  zero  if  this  is  unknown; 

uxyt_x  returns  real  value  of  ux(a;,y,t)  or  zero  if  this  is  unknown; 

uxyt_y  returns  real  value  of  uy(x,y,t)  or  zero  if  this  is  unknown; 

vxyt  returns  real  value  of  vertical  displacement  v(x^y}t)  or  zero  if  this  is  unknown; 

vxyt_x  returns  real  value  of  vx(x,y,t)  or  zero  if  this  is  unknown; 

vxyt.y  returns  real  value  of  %(x,y,£)  or  zero  if  this  is  unknown; 

f  lxyt  returns  real  value  of  horizontal  body  force  /i(^,y,  t); 


92 


CHAPTER  9 .  OBTAINING  AND  USING  THE  SOFTWARE 


f2xyt  returns  real  value  of  vertical  body  force  /2(:r,y,£); 

Sigmallxyt  returns  real  value  of  direct  stress  <7n(:r,y,  t)  or  zero  if  unknown; 
Sigma22xyt  returns  real  value  of  direct  stress  a22(z,y,  t)  or  zero  if  unknown; 
Sigmal2xyt  returns  real  value  of  shear  stress  cr12(a;,y,t)  or  zero  if  unknown; 
Epsilonllxyt  returns  real  value  of  direct  strain  £\\{x,y,t)  or  zero  if  unknown; 
Epsilon22xyt  returns  real  value  of  direct  strain  S22(x,y,t)  or  zero  if  unknown; 
Epsilonl2xyt  returns  real  value  of  shear  strain  £12 (x>y,t)  or  zero  if  unknown; 
glxyt  returns  real  value  of  horizontal  traction  yi(a;,y,£); 
g2xyt  returns  real  value  of  vertical  traction  y2(:r,  y,t) ; 

The  other  functions  in  the  source  code  are  connected  with  the  relaxation  functions  A  and 
ji  and  can  be  accepted  as  they  stand. 

In  our  trestle  example  the  only  non-zero  traction  is  vertical,  of  of  —lMPa  in  mag¬ 
nitude,  and  applied  along  only  one  edge.  Thus  we  make  sure  that  the  function  glxyt 
contains  the  lines, 

#elif  defined(TRESTLE) 
total  =  (real) 0.0; 

#else 


since  there  are  no  horizontal  tractions,  and  also  that  the  function  g2xyt  contains 


#elif  defined (TRESTLE) 
total  =  -(real)1.0e6; 
#else 


(The  type  cast  (real)  is  not  mandatory  in  these  statements — it  can  be  ignored  and  left 
out  if  desired.)  Notice  that  we  do  not  have  to  specify  that  y2  is  only  non-zero  for  the 
specific  coordinates  (z,y)  occurring  on  edge  0 — this  was  taken  care  of  with  the  boundary 
condition  codes  supplied  to  the  mesh  generator.  Since  this  is  the  only  edge  with  the  non¬ 
zero  traction  boundary  code  3,  this  C  function  is  only  invoked  for  element  edges  occurring 
on  this  boundary  edge. 

The  remaining  task  is  to  define  the  material  properties.  At  the  moment  the  code  is  still 
set  to  use  the  values  given  by  (7.13)  to  (7.16).  However,  to  use  the  more  realistic  values  for 
Mar  any  1  given  in  Section  8.5  we  need  to  make  a  small  change  to  another  #def  ine  directive. 
This  one  is  near  the  top  of  the  source  file  seedproj/quasivis/native/source/viscous  .  c. 
Edit  this  file  and  change  the  line 
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#def ine  EXACT 


to 

#def ine  MARANYL 


Upon  compilation,  the  declarations  at  the  head  of  this  code  now  set  up  two  structure-type 
variables  that  describe  the  Maranyl  Nylon  6.6  relaxation  function. 

We  are  now  ready  to  run  the  code.  Change  to  the  seedproj/quasivis/data directory 
and  list  the  files.  You  should  see  the  directories  exact,  crack,  leverarm,  web  and  our 
new  one  trestle,  and  also  two  /sbin/sh  shell  scripts  adapt  and  run.  Do  not  attempt 
to  use  the  latter,  it  was  written  (rather  badly)  in  order  to  automate  the  generation  of  the 
numerical  results  given  earlier.  Instead  we  shall  use  the  simpler  adapt  script.  Type, 

adapt  trestle  1.5 

This  moves  all  working  files  into  the  dump  directory;  sets  the  $PATH  variable  so  the  shell 
can  find  the  executables  advance,  Xwin  and  fem;  generates  the  mesh  according  to  the 
* .  inf  and  * .  dat  files  in  the  trestle  directory;  makes  and  links  any  updated  source  code; 
and,  finally  launches  the  finite  element  solver  fem  with  the  command, 

fem  -X  -memscale  100  -adapt  $T0L  1.0 

Here  the  argument  -X  requests  the  X- window  graphics,  while  -memscale  100  allocates 
memory  for  a  mesh  a  hundred  times  more  dense  that  the  starting  mesh  (the  one  given 
earlier  on  the  left  of  Figure  9.3).  The  argument  -adapt  tells  the  code  to  solve  in  the 
adaptive  mode  with  tolerance  TOL  given  by  the  third  command  line  argument  to  the 
adapt  script.  In  this  case  the  environment  variable  $T0L  has  the  value  TOL  =  1.5.  The 
final  value  1.0  tells  the  code  to  apply  100%  of  this  tolerance  to  spatial  error  control  and 
none  to  temporal  error  control.  Eventually,  when  some  form  of  temporal  adaptivity  is 
implemented,  it  will  be  possible  to  have  a  command-line  option  to  fem  taking  the  form, 

-adapt  1.2  0.4 

This  would  specify  a  global  tolerance  of  TOL  =  1.2  with  only  40%  (i.e.  TOLq  =  0.48)  going 
toward  controlling  the  spatial  errors  via  the  residual  term  £n-  However,  until  temporal 
adaptivity  is  implemented  the  second  value  for  the  -adapt  flag  should  always  be  1.0. 

Executing  the  command 

adapt  trestle  1.5 
should  produce  the  output  lines 

ki  N  *l«a  I  lu  I  |_E  ||uh||_E  |uhM«l  |uh|+«st|a|  f|u-uhj|_E  lit  ||«||_E  inf|u-uh|  Mu-uhll  ||»-*h|| 

1.0000*+00  1123  0.000a+00  1.936a+01  2.738a+01  1.939a+01  1.935795a+0l  1.138819*400  3.343«-03  9.364a-06'  2.670*405 
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Figure  9.4:  Adapted  meshes  for  the  trestle  problem  with  TOL^  =  1.5  on  the  left  and 
TOLft  =  1.0  on  the  right. 


(although  you’ll  need  a  wide  terminal  to  see  this  clearly),  as  well  as  the  adapted  mesh 
shown  on  the  left  of  Figure  9.4. 

To  see  the  effect  of  the  adaptivity  reduce  the  tolerance  to  1.0  and  execute  again  with 
the  line 

adapt  trestle  1.0 
This  produces  the  output  lines, 

ki  N  «l«m  I  |ul  I  _E  ||  uh  1 1  _E  |uh|  +  |*|  |uh|4«rt|a|  ||u-uh||_E  IUILE  influ-uhl  llu-uhll  1 1  ®-«h  |  | 

1.0000«+00  2097  0.000«+00  1.955*401  2.764*401  1.957*401  1.954769*401  8.495796«~01  3.41i*-03  9.558«-05  2.695*405 

and  the  adapted  mesh  shown  on  the  right  of  Figure  9.4. 

All  of  the  program  execution  takes  place  in  the  . .  /dump/  directory,  and  here  you 
will  find  four  files  of  interest.  The  first  is  mesh0.tex  which  is  an  example  of  the  input 
files  used  to  illustrate  the  displaced  meshes  in  this  document.  It  is  a  mixture  of  WFftK.  2e 
source  which  uses  the  TjgXdraw  graphics  package.  The  other  three  files  are  sigmall_0.m, 
sigma22_0.m  and  sigmal2_0.m.  These  are  Matlab  script  files  which  may  be  executed 
within  Matlab  as  they  stand  and  will  produce  the  surface  plots  of  the  stresses.  These  were 
also  used  earlier  in  this  document  to  illustrate  the  numerical  solutions. 

Note  that  the  directory  dump  is  named  deliberately.  All  files  in  it  are  disposable,  in 
that  they  can  regenerated,  and  are  potentially  large  and  so  can  eat  up  disk  space.  It  is 
recommended  that  this  directory  be  cleaned  out  once  you  have  finished  using  the  code. 
However,  leave  the  (empty)  dump  directory  in  place  since  its  existence  is  assumed  and 
required  by  the  script  files  adapt  and  run. 


9.6  Some  comments 

There  is  clearly  room  for  a  great  deal  of  improvement  in  the  way  the  code  is  currently 
constructed  and  configured.  A  possible  way  forward  here  would  be  to  provide  a  platform- 
independent  Java  GUI  which  could  invoke  the  existing  native  C  code  through  a  native 
Java  method  of  a  larger  object.  This  is  however  not  research  within  the  scope  of  the  seed 
project  and  so  we  suggest  it  only  as  a  possibility  for  the  future. 


Chapter  10 


Suggestions  for  further  work 


10.1  Overview 

In  this  short  concluding  chapter  we  give  brief  details  of  the  code  development  work  carried 
out  between  ourselves  and  Dr.  A.R.  Johnson  during  his  visit  to  BICOM  in  March  1998. 
We  then  give  a  list  of  research  areas  that  could  follow  naturally  from  the  work  detailed  in 
this  report. 


10.2  A.R.  Johnson’s  visit  to  BICOM 

During  the  first  week  of  March  1998  Dr.  A.R.  Johnson  (ARL,  VTC,  NASA,  Langley,  VA, 
USA)  visited  BICOM  with  the  aim  of  investigating  how  adaptivity  could  be  built  in  to  his 
existing  work  on  internal  variables  (see  for  example  [24,  25]).  During  this  short  time  we 
and  Dr.  Johnson  were  able  to  generate  a  Fortran  FUNCTION  that  was  easily  added  on  to  one 
of  his  existing  finite  element  codes,  in  this  case  for  solving  viscoelastic  beam  problems. 
This  FUNCTION  adaptively  and  automatically  solves  the  internal  variable  equation  over 
each  of  the  time  steps  taken  by  the  larger  finite  element  solver.  Its  modular  nature  means 
that  it  can  also  be  easily  transplanted  to  function  with  codes  written  to  solve  problems 
involving  structures  other  than  beams.  In  this  section  we  will  briefly  outline  the  nature 
of  the  algorithm  and  the  type  of  error  control  produced,  and  also  give  a  “listing”  of  the 
FUNCTION  to  illustrate  how  this  process  can  be  easily  encapsulated  in  a  simple  module. 

In  [24,  Equation  20.38]  Johnson  and  Tessler  show  how  a  finite  element  discretization 
of  a  quasistatic  viscoelasticity  problem  can  be  written  in  the  form, 

Ku  =  F mech  -F1 vise  at  time  t.  (10. 1) 

Here:  u  is  the  vector  of  unknown  nodal  displacements;  K  is  the  finite  element  stiffness 
matrix;  Fmech  the  vector  of  external  loads;  and,  FViSc  a  vector  of  viscous  forces  containing 
the  viscoelastic  effects.  The  vector  FviSC  is  obtained  by  assembling  the  local  element  level 
viscous  forces  via, 

F  visC  =  ^  k  li, 
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where  *k  is  an  element-level  viscous-stiffness  matrix,  and  *u  a  vector  of  unknown  internal 
displacements  defined  at  the  element  level. 

As  shown  by  Johnson  and  Tessler  in  [24,  Equation  20.37]  these  internal  displacements 
are  given  on  each  element  by  the  solution  of  the  evolution  equation, 


d  *u  *u 
dt  ^  r 


du 
dt  ’ 


(10.2) 


Thus,  to  find  the  displacement  u(t)  at  the  given  time  t  one  has  only  to  solve  (10.1)  to 
give, 


u(t)  —  K  (^F  mech  -^Visc^* 

It  is  natural  to  assume  that  -Fmech  is  known  precisely  and  so  the  time  discretization  error 
in  this  equation  is  effectively  present  in  the  term  Fv jsc.  It  is  therefore  important  to  have 
an  adaptive  solution  algorithm  for  the  evolution  equation  (10.2)  so  that  the  viscous  loads 
can  be  built  accurately  and  cheaply.  It  was  the  derivation  of  such  an  algorithm,  and  its 
realization  and  modular  implementation  as  a  Fortran  FUNCTION,  that  formed  the  focus  of 
the  collaborative  work  carried  out  at  BICOM  during  Dr.  Johnson’s  visit. 

We  give  brief  details  of  this  work  below,  and  we  finish  be  suggesting  ways  in  which 
this  work  could  be  extended. 

The  internal  variable  equation  can  be  considered  as  a  collection  of  scalar  ODE  prob¬ 
lems,  each  having  the  form:  find  u  such  that, 

^  -  in  [0,T],  with  u(0)  =  uo,  (10.3) 

dt  r  r 

and  where  T  >  0,  u0,  g  and  r  >  0  are  given  data.  Here,  of  course,  [0,T]  represents  a 
single  time  interval  within  the  time  stepping  loop  of  the  overall  finite  element  solver,  and 
g  is  derived  from  the  global  displacement  vector  u.  Below  we  set, 

g  _  du 
r  dt  ’ 

although  we  could  also  work  in  terms  of  reversed  internal  variables  with  essentially  no 
change  to  the  adaptive  algorithm. 

Our  approach  to  deriving  a  numerical  scheme  is  based  on  discretizing  the  single  time 
interval  [0,  T]  into  subintervals  in  the  usual  way, 

0  =  t0  <  t\  <  •  •  •  <  tq  <  •  •  •  <  tN  =  r, 

and  then  defining  the  time  steps  kq  =  tq  -  tq- i,  and  subintervals  Jq  (tq- 1,  tq).  Before 
proceeding  we  note  that  it  is  possible  to  “solve”  (10.3)  explicitly  by  introducing  an  inte¬ 
grating  factor,  but  one  would  then  need  to  evaluate  the  resulting  integral  with  “sufficient” 
accuracy.  The  approach  described  below  could  thus  also  be  described  as  an  adaptive 
quadrature  algorithm,  and  in  particular  can  still  be  applied  with  essentially  no  changes  in 
the  nonlinear  case  where  r  =  r{u). 
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Our  discretization  of  the  internal  variable  equation  (10.3)  is  based  on  the  finite  ele¬ 
ment  method.  We  introduce  a  piecewise  constant  (i.e.  constant  on  each  Jq)  approxima¬ 
tion  U  to  u  and,  in  particular,  denote  by  Uq  the  constant  approximation  to  u  during  times 
t  e  Jq.  Then,  the  Galerkin  finite  element  approximation  of  (10.3)  can  be  written  as, 

(1  +  kq/r)Uq  =  Uq-i  +  f  q  —  dt  with  U0  =  u0 . 

Jtq- 1  T 

One  obtains  this  time  stepping  scheme  by  taking  the  scalar  product  of  (10.3)  with  a 
piecewise  constant  test  function,  replacing  u  with  U  (interpreting  the  derivative  in  a 
distributional  sense)  and  then  applying  standard  finite  element  methodology. 

Using  duality  arguments  we  then  arrive  at  the  a  posteriori  error  estimate, 

\u(tN)  -  UN I  <  (1  -  e-Wr)  ma  (~||5  -  U\\LM  +  I  Uq  . 

Note  that  the  quantities  on  the  right  are  completely  determined  by  the  given  data  and  the 
computed  solution,  and  so  this  bound  can  form  the  basis  of  an  adaptive  algorithm.  We 
describe  this  more  fully  below. 

The  drawback  with  the  approach  just  taken  is  that  the  piecewise  constant  approxima¬ 
tion  to  u  is  not  particularly  accurate.  To  address  this  we  also  developed  and  implemented 
a  continuous  piecewise  linear  finite  element  approximation  to  (10.3),  where  we  retained 
the  piecewise  constant  test  function.  In  this  case  the  time  stepping  scheme  reduces  to, 

(1  +  kq/2T)Uq  =  {1-  kq/2T)Uq—\  +  [ k  dt. 

Jtq-l  T 

For  this  scheme  we  derived  the  a  posteriori  error  estimate, 

N<jv)  -  <  (1  -  e~tN/T)  maxN  {^||r||Loo(Js)}, 


where, 


is  the  residual  and  is  computable. 

We  use  this  a  posteriori  error  estimate  to  generate  an  adaptive  time  stepping  scheme 
in  the  following  way.  Suppose  we  want  to  compute  a  numerical  solution  for  which, 

M*N)-U(tN)\  <TOL 

is  guaranteed,  where  TOL  >  0  is  a  user-specified  tolerance  level,  then  it  is  sufficient  to 
ensure  that, 


(1  -  e  tN/T)  max^  {fc9||r-||LooW)}  =  TOL. 


98 


CHAPTER  10.  SUGGESTIONS  FOR  FURTHER  WORK 


This  in  turn  is  also  guaranteed  if  we  ensure  that, 

MHIw^TOLa-e-^)-1 


for  every  q.  Rearranging  this  leads  to  a  simple  rule  for  selecting  the  time  steps  in  an 
iterative,  or  adaptive,  manner: 


£  new  __ 


T0L(1  -  e~tj'f/'T)~1 

IMUooW) 


This  adaptive  solver  for  (10.3)  can  be  encapsulated  as  a  Fortran  FUNCTION  as  listed  below. 


C********************************************************* 

c  **********  CONTINUOUS  ********** 

c  **********  PIECEWISE  LINEAR  APPROX  ********** 

c -  MUST  be  declared  as  REAL*8  - 

c  usage: 

c  ustar  =  ivadaptl (du/dt ,  ic,  tau.n,  t_{m-l},  t_m,  TOL,  steps,  mk) 
c  minimum  predicted  time  step  may  be  returned  in  mk  if  mk  is  larger 
c  on  entry.  Otherwise  the  value  in  mk  is  not  altered 
c  Note  that  "steps"  return  the  number  of  steps  taken  in  the 
c  time  interval  (tlo,  thi) . 
c 

real*8  function  ivadaptl (rhs,ic,tau, tlo, thi, TOL, steps, mk) 


c 


c 


c 


c 

c 


real*8  rhs,  ic,  tau,  tlo,  thi,  TOL,  t,  mk 

time  step,  new  time  step,  previous  and  current  iterates 

real*8  k,  new.k,  prevu,  curru 

real*8  numer,  denom  !  for  time  step  control 

real*8  denoml,  denom2  !  for  time  step  control 

integer  iter ,  steps 


steps  =  0 
prevu  =  ic 
k  =  thi  -  tlo 
t  =  tlo 


!  record  the  number  of  steps  taken 
!  set  initial  condition 
!  guess  initial  time  step 
!  current  time  level  is  (t,  t+k) 


NOTE:  next  quantity  can  be  supplied  to  save  "exp"  calculations 
numer  =  TOL  /  (  l.OdO  -  exp(- (thi-tlo) /tau)  )  !  for  adaptivity 
do  iter  =  1,  1000  !  no  idea  how  many  steps  are  needed! 

steps  =  steps  +  1 
get  next  solution 

curru  =  (1  -  k/2 . OdO/ tau) *prevu  +  rhs*k 

curru  =  curru  /  (1  +  k/2.0d0/tau) 

test  error  condition  by  computing  a  new  time  step, 

first  find  the  max.  value  in  the  denominator 


denoml  =  abs(  rhs  -  prevu/tau  -  (curru-prevu)/k  ) 

denom2  =  abs(  rhs  -  curru/tau  -  (curru-prevu) /k  ) 

denom  =  max (denoml,  denom2)  !  next  line  can  throw  an  exception 

new_k  =  numer  /  denom  !  predict  a  new  time  step 

!  and  act  accordingly 

!  remember  minimum  predicted  time  step 


mk  =  min(mk,  new_k) 
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if(  new_k  .GE.  k  )  then  !  advance  to  next  time  level 

t  =  t  +  k 

k  =  MIN(  new_k,  thi  -  t  ) 
prevn  =  curru 

if(  t  .GE.  thi  )  goto  1  !  break  out  if  we  have  reached  final  time 


else  !  recompute  this  time  level 

k  =  new_k  !  reduce  k  and  try  again  at  this  level 

end  if 
enddo 

write(*,*)  3 . ERROR  in  ivadaptlQ  . , 


1  ivadaptl  =  curru  !  return  solution  at  thi 
end 

C***  *********  j(t^*jfc*jf:^****J|cj|e3jcjJc)Jcj|c*j(cJ(c***3|c*************J»c***************4:****J»c**>tc 


The  initial  condition  is  given  by  ic  and  the  right  hand  side  g/r  is  supplied  as  the 
constant  (on  the  time  step)  rhs.  The  interval  [0,T]  is  supplied  through  the  argument  list 
in  the  more  general  form  and  the  FUNCTION  returns  the  approximation  U(tm) 

to  u(tm). 

This  approach  could  easily  be  modified  to  develop  higher  order  solvers  for  the  internal 
variable  equation,  and  then  applied  to  more  or  less  arbitrary  viscoelasticity  problems  that 
can  be  formulated  with  such  evolution  equations.  The  approach  is  also  highly  suited  to 
constitutively  nonlinear  problems  modelled  with  a  reduced  time.  In  such  a  case  r  =  r(u), 
but  this  will  introduce  no  essential  complications  into  the  adaptive  scheme. 


10.3  Additional  research  areas 

This  section  is  deliberately  short  and  “bullet  pointed”.  We  want  only  to  suggest  areas 
that  could  be  fruitfully  investigated  on  the  back  of  this  work,  and  not  pre-empt  the  fine 
detail  of  the  topics. 

•  In  a  time  dependent  problem  it  is  crucial  for  an  adaptive  algorithm  to  be  able  to 
selectively  de-refine  as  well  as  refine  the  space  mesh  as  features  in  the  solution  evolve 
or  decay  over  time.  As  we  explained  earlier  in  Chapter  8,  our  current  a  posteriori 
error  estimates  (given  in  detail  in  [55,  56])  do  not  allow  de-refinement  without  the  in¬ 
clusion  of  a  complicated  and  expensive  residual  term.  We  believe  the  internal  variable 
algorithm  developed  in  this  report  will  allow  remove  this  expensive  term  for  good  and 
therefore  allow  mesh  de-refinement  with  only  a  modest  amount  of  additional  expense. 
We  plan  to  investigate  this  further  at  a  later  date. 

•  Temporal  error  control  should  be  achievable  by  measuring  the  energy  error  in  an 
appropriately  weak  norm.  Results  for  a  prototype  problem  have  already  been  given  in 
[53],  and  the  extension  of  this  work  to  the  space-time  quasistatic  problem  is  underway 
in  [56]. 

•  Ultimately  the  goal  must  be  solve  physically  realistic  nonlinear  problems.  We  feel 
that  extending  our  results  to  constitutively  nonlinear  problems,  characterized  by  a 
reduced  time,  would  be  fairly  straightforward,  and  it  may  even  be  true  that  the 
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stability  estimates  in  [57]  will  hold  without  modification.  However,  nonlinearity  aris¬ 
ing  from  large  deformations  is  a  much  more  challenging  problem  and  as  yet  we  have 
no  theoretical  insights  that  will  help  in  the  construction  of  adaptive  software. 

•  Quasistatic  linear  plate  or  shell  problems  could  also  be  addressed  in  essentially  the 
same  way  as  described  in  this  report  for  the  “standard”  two-  and  three-dimensional 
problem.  In  particular,  the  stability  estimates  in  [57]  should  apply  directly. 

•  The  space-time  finite  element  method  as  illustrated  in  this  report  can  also  be  used  to 
generate  a  posteriori  error  estimates  and  adaptive  algorithms  for  dynamic  viscoelas¬ 
ticity  problems.  We  hope  to  begin  to  look  at  this  subject  soon. 


Part  IV 
References 
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